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thermophysiea!  problems  and  prepare  them  for  program  input.  In  addition,  the  program  contains 
numerous  subroutines  for  handling  interrelated  complex  phenomena  such  as  sublimation;  diffuse 
radiation  within  enclosures;  simultaneous,  one-dimensional,  incompressible  fluid  flow  including 
valving  and  transport  delay  effects;  etc.  It  can  handle  all  three  types  of  heat  transfer— conduction, 
convection,  and  radiation.  The  optional  combinations  of  these  capabilities,  in  conjunction  with 
the  model  size  allowable  (4000  nodes  on  a  65k-core  machine),  make  CINDA  an  extremely  potent 
analytical  tool  for  thermal  systems  analysis  in  the  hands  of  a  competent  engineer  analyst.  Its  uses 
include  determining  temperatures  of  structures  such  as  bridges,  rockets,  and  buildings;  finding 
pooling  requirements  for  electric  circuits;  and  studying  the  thermal  properties  of  adverse  thermal 
systems  such  as  nuclear  reactors  and  automobile  engines. 

The  programs  on  pp.  104, 105, 106  are  adaptations  of  similar  programs  published  in  “MITAS 
User  Information  Manual  (Martin  Marietta  Thermal  Analyzer  System),”  CYBERNET 
Publications  Division.  Control  Data  Corporation,  Copyright  1972. 
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I.  INTRODUCTION 
Background 

When  it  was  recognized  that  a  thermal  analyzing  system  was  needed  at  the  Naval 
Research  Laboratory,  the  Chrysler  Numerical  Differencing  Analyzer  for  third-generation 
computers  (CINDA-3G)*  was  obtained  for  use  on  the  CDC-3800  computer.  CINDA-3G 
was  developed  by  the  Thermodynamic  Section  of  the  Aerospace  Physics  Branch  of  the 
Chrysler  Corporation  Space  Division  at  NASA’s  Michoud  Assembly  Facility.  A  major 
portion  of  this  work  was  done  as  a  NASA-funded  project  from  the  Manned  Spacecraft 
Center  in  Houston,  Texas. 

CINDA  was  programmed  in  FORTRAN  V  for  the  Univac  1108.  Because  of  differ¬ 
ences  between  the  computers,  modifications  to  the  program  were  necessary  before  it 
could  be  usable  on  the  CDC-3800.  The  task  of  conversion  required  compensating  for 
compiler  differences  (FORTRAN  V  vs  3800  FORTRAN),  rewriting  some  routines  in 
COMPASS,  and  in  general,  adapting  CINDA  to  the  Drum  SCOPE  system.  Documenta¬ 
tion  supplied  v/ith  the  Univac  version  of  CINDA  was  a  most  useful  aid  in  the  conversion 
process.  This  NRL  report  represents  the  CDC-3800  version  of  CINDA-3G,  and  is  partly 
a  rewrite  of  the  original  Chrysler  document.  While  most  of  the  sections  have  been  only 
slightly  modified,  the  section  on  tape  usage  and  deck  setups  is  strictly  applicable  to 
CDC-3800  software.  The  external  plotting  package  in  Section  VI  for  the  CALCOMP 
plotter  replaces  plotting  routines  used  by  the  Univac  1108  for  the  SC-4020. 


Overview 

This  programming  manual  deals  with  the  CINDA  system  in  two  general  categories— 
the  logic  and  data  needed  in  setting  up  the  problem  data  deck,  and  the  actual  structure 
of  the  operating  system. 

The  logic  in  constructing  a  problem  for  CINDA  involves  developing  a  lumped- 
parameter  representation  of  the  physical  problem.  This  model  simulates  the  elements  of 
heat  transfer,  and  the  user  must  supply  the  corresponding  network  data,  which  will  be 
used  by  one  or  more  routines  selected  from  a  large  subroutine  library.  The  user  must 
determine  which  routines  are  needed  and  the  order  in  which  they  are  to  be  activated. 
This  information  and  the  other  related  logic  are  entered  in  a  modified  FORTRAN 
language.  The  major  routines  involved  use  various  iterative  techniques  for  solutions, 
with  the  program-formed  compute  sequences  minimizing  the  required  operations. 


Note:  Manuscript  submitted  September  5,  1973. 

♦D.R.  Lewis,  J.D.  Gaski,  and  L.R.  Thompson,  “Chrysler  Improved  Numerical  Differencing  Analyzer  for 
3rd  Generation  Computers,"  Technical  Note  TN-AP-G7-287,  Chrysler  Corporation  Space  Division,  New 
Orleans,  La.,  1967. 
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CINDA  is  not  merely  a  single  execution  program,  but  is  an  operating  system  in  itself. 
It  consists  of  two  software  packages  of  its  own  (the  preprocessor  and  subroutine  library 
file)  and  is  also  quite  dependent  on  the  computer’s  system  software  (at  NRL,  the  CDC- 
3800’s  Drum  SCOPE  system).  This  dependence  is  largely  due  to  language  and  allocating 
differences  of  FORTRAN  compilers,  as  well  as  computer  work  size  and  assembly  language 
routines.  The  preprocessor  reads  and  processes  the  two  types  of  user  data— the  network 
data  and  the  logic  data.  From  the  former,  its  output  is  a  binary  data  file,  and  from  the 
latter,  it  generates  five  FORTRAN  subroutines  which,  when  compiled,  are  referred  to  as 
the  processor.  The  processor  makes  user-requested  calls  to  subroutines  in  the  library  file, 
executing  the  binary  data  output  from  the  preprocessor.  The  user’s  logic  controls  the 
processor  until  the  end  of  the  job.  For  any  particular  problem  the  arrays  in  the  processor 
are  dimensioned  exactly  as  needed  by  the  preprocessor.  This  feature,  together  with  user- 
controlled  output,  saves  computer  time  and  money. 


II.  DISCUSSION 
Lumped-Parameter  Representation 

The  key  to  utilizing  a  network-type  analysis  program  lies  in  the  users’  ability  to  de¬ 
velop  a  lumped-parameter  representation  of  the  physical  problem.  Once  this  is  done, 
superimposing  the  network  mesh  is  a  mechanical  task  at  most  and  the  numbering  of  the 
network  element.,  is  simple  although  perhaps  tedious.  It  might  be  said  that  the  network 
representation  is  a  “crutch”  for  the  engineer,  but  it  does  simplify  the  data  logistics  and 
allow  easy  preparation  of  data  input  to  the  program.  In  addition,  it  allows  the  user  to 
identify  uniquely  any  element  in  the  network  and  modify  its  value  or  function  during 
the  analysis  as  well  as  sense  any  potential  or  current  flow  in  the  network.  Another  fea¬ 
ture  of  the  network  is  that  it  has  a  one-to-one  correspondence  to  the  mathematical  model 
as  well  as  the  physical  model. 

Perhaps  the  most  critical  aspect  of  the  lumped-parameter  approach  is  determining 
the  lump  size.  There  are  methods  for  optimizing  the  lump  size,  but  they  usually  involve 
more  analytical  effort  and  computer  time  than  the  original  analysis.  One  must  also  keep 
in  mind  that  for  a  transient  problem,  time  is  being  lumped  as  well  as  space.  Of  prime 
importance  is  what  information  is  being  sought  from  the  analysis.  If  spot  temperatures 
are  being  sought,  nodes  must  at  least  fall  on  the  spots  and  not  include  much  more  physi¬ 
cally  than  would  be  expected  to  exist  at  a  relatively  similar  temperature.  Nodes  must 
fall  at  end  points  when  a  temperature  gradient  is  sought.  Of  necessity,  lumping  must  be 
fairly  fine  where  isotherms  are  sought.  Lumping  should  be  coarse  in  areas  of  high  ther¬ 
mal  conductivity.  When  nonlinear  properties  are  being  evaluated,  the  lumping  should  be 
fine  enough  so  that  extreme  gradients  are  not  encountered.  The  lumping  is  also  depend¬ 
ent  on  the  severity  of  the  nonlinearity. 

To  reduce  round-off  error,  the  explicit  stability  criteria  of  the  lump  (the  capacitance 
value  divided  by  the  summation  of  conductor  values  into  the  node)  should  be  held  fairly 
constant.  The  value  C/XG  is  directly  proportional  to  the  square  of  the  distance  between 
nodes.  Although  refining  the  lumped-parameter  representation  will  yield  more  accurate 
answers,  halving  the  distance  between  nodes  decreases  the  stability  criteria  by  a  factor  of 
four  and  increases  the  numner  of  nodes  by  a  factor  of  two,  four,  or  eight  depending 
upon  whether  the  problem  is  one,  two,  or  three  dimensional.  For  the  explicit  case, 


2 


NP.L  REPORT  7656 


halving  the  distance  between  nodes  increases  the  machine  time  for  transient  analysis  by  a 
factor  of  8,  16,  or  32,  respectively.  The  increase  in  solution  time  for  the  implicit  methods 
is  somewhat  less  but  proportional. 

When  lumping  the  time  space,  consideration  must  be  given  to  the  frequency  of  the 
boundary  conditions.  A  time  step  must  not  step  over  boundary  excitation  points  or 
they  will  be  missed.  Do  not  step  over  pulses;  rather,  rise  and  fall  with  them.  Generally 
the  computation  interval  for  the  explicit,  methods  is  sufficiently  small  so  that  frequency 
effects  can  be  ignored.  However,  care  must  be  exercised  when  specifying  the  time  step 
for  implicit  methods.  If  only  a  small  portion  of  a  transient  analysis  involves  frequency 
considerations,  the  time  step  used  may  be  selectively  restricted  for  that  interval.  By 
setting  the  maximum  time  step  allowed  as  a  function  of  time,  we  may  utilize  an  inter¬ 
polation  call  to  vary  it  accordingly. 

One  must  also  realize  that  the  problem  being  solved  is  linearized  over  the  time  step. 
Heating  rate  calculations  are  usually  computed  for  a  time  point  and  then  applied  to  a 
time  space.  If  the  rates  are  nonlinear,  a  certain  amount  of  error  is  introduced,  partic¬ 
ularly  so  with  radiation.  These  nonlinear  effects  may  cause  almost  any  method  of  solu¬ 
tion  to  diverge.  A  brute  force  method  for  forcing  convergence  is  to  limit  the  temperature 
change  allowed  over  the  time  space.  Consideration  of  the  factors  mentioned  above,  cou¬ 
pled  with  some  experience  in  using  the  program,  will  aid  the  observant  analyst  in  choosing 
lump  sizes  that  will  yield  answers  of  sufficient  engineering  accuracy  with  a  reasonable 
amount  of  computer  time. 

Figure  1  shows  the  lumped-parameter  representation  and  network  superposition  of  a 
one-dimensional  heat  transfer  problem. 


Fig.  1— One-dimensional  network 


The  “node”  points  are  centered  within  the  lumps,  and  temperatures  at  the  nodes  are  con¬ 
sidered  uniform  throughout  the  lump.  The  capacitors  hung  from  the  nodes  indicate  the 
ability  of  the  lump  to  store  thermal  energy.  Capacitance  values  are  calculated  as  lump 
volume  times  density  times  specific  heat.  The  conductors  (electrical  symbol  G)  represent 
the  capability  for  transmitting  thermal  energy  from  one  lump  to  another.  Conductor 
values  for  energy  transmission  through  solids  are  calculated  as  thermal  conductivity  times 
the  energy  cross-sectional  flow  area  divided  by  path  length  (distance  between  nodes). 
Conductor  values  for  convective  heat  transfer  are  calculated  as  the  convection  coefficient 
times  the  energy  cross-sectional  flow  area.  Conductors  representing  energy  transfer  by 
radiation  are  usually  indicated  by  crossed  arrows  over  the  conductor  symbol.  Radiation 
transfer  is  nonlinear;  it  is  proportional  to  the  difference  of  the  absolute  temperatures 
raised  to  the  fourth  power.  Utilization  of  the  Farenheit  system  allows  easy  automation 
of  this  nonlinear  transfer  function  by  the  program  and  reduces  the  radiation  conductor 
value  to  the  product  of  the  Stephan-Boltzmann  constant  times  the  surface  area  times  the 
net  radiant  interchange  factor  (script  F). 
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Basics  of  Finite  Differencing 


The  oncept  of  network  superposition  on  the  lumped-parameter  representation  of  a 
physical  system  is  easy  to  grasp.  Describing  the  network  to  the  program  is  also  quite 
straightforward.  Having  described  a  network  to  the  program,  what  information  have  we 
really  supplied  and  what  does  the  program  do  with  it?  Basically,  we  desire  the  solution 
to  a  simultaneous  set  of  partial  differential  equations  of  the  diffusion  type;  i.e., 


~  =  «V2T  + 
<)t 


S, 


ii 

0y2 


(1) 


That  the  diffusivity  (a  =  k/pCp)  may  l)e  temperature  varying  or  nonlinear  radiation  trans¬ 
fer  occurring  is  immaterial  at  this  point.  Of  importance  is  how  Eq.  (1)  is  finite  differ¬ 
enced  and  its  relationship  to  the  network  and  energy  flow  equations  more  commonly 
utilized  by  the  engineer.  The  partial  of  the  T-state  variable  with  respect  to  time  is  finite 
differenced  across  the  time  space  as  follows: 


3T  _  T'  -  T 
3t  At  ’ 


(2) 


where  the  prime  indicates  the  new  T  value  after  passage  of  the  At  time  step. 

The  right  side  of  Eq.  (1)  could  be  written  with  T  primed  to  indicate  implicit  “back¬ 
ward”  differencing  or  unprimed  to  indicate  explicit  “forward”  differencing.  The  follow¬ 
ing  equation  is  illustrative  of  how  “backward”  and  “forward”  combinations  may  be 
obtained. 


^  =  0(aV2T  +  S)  +  (1  -  0)(a'V2T'  +  S') ,  0  <  p  <  1  .  (3) 

at 

Any  value  of  P  less  than  1  yields  an  implicit  set  of  equations  which  must  be  solved  in  a 
simultaneous  manner  (more  than  one  unknown  exists  in  each  equation).  Any  value  of  p 
equal  to  or  less  than  1/2  yields  an  unconditionally  stable  set  of  equations  or  in  other 
words,  any  time  step  desired  may  be  used.  Values  of  p  greater  than  1/2  invoke  stability 
criteria  or  limitations  on  the  magnitude  of  the  time  step.  A  value  of  P  equal  to  1/2 
yields  an  unconditionally  stable  implicit  set  of  equations  commonly  known  as  “forward- 
backward”  differencing  or  the  Crank-Nicholson  method.  Various  transformations  or  first 
order  integration  applied  to  Eq.  (1)  generally  yield  an  implicit  set  of  equations  similar  to 
Eq.  (3)  with  p  equal  to  1/2.  The  following  finite  difference  approach  generally  applies 
to  transformed  equations. 

Let  us  consider  the  right  side  of  Eq.  (3)  with  p  -  1  and  rewrite  it  as  follows: 


arV2T  +  S 


a_  fIL  _  JT)  +  /il_  _  i!L\  +  -«  fdT 

Ax  \3x-  3x+/  Ay  \3y-  3y  +/  Az  \3z- 


+  S 


(4) 


The  minus  or  plus  signs  on  the  first  partial  terms  indicate  that  they  are  taker,  on  the 
negative  or  positive  side,  respectively,  of  the  point  under  consideration  and  always  in  the 
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same  direction.  If  we  consider  three  consecutive  points  (1,  2,  and  3)  ascending  in  the  x 
direction  we  can  complete  the  finite  difference  of  the  x  portion  of  Eq.  (4)  as  follows: 


a 

Ax 


(W* 

_  ^ 

a  /T,  -  T2  ]  T3  -  T2\ 

\3x- 

0X+/ 

Ax  \  Ax-  Ax+  / 

'5) 


Applying  the  above  step  to  the  y  and  z  portions  of  Eq.  (4)  yields  the  common  denomi¬ 
nator  of  volume  (V  =  Ax*  Ay*  Az).  Using  Eq.  (3)  with  0=1,  finite  differencing  with 
the  steps  used  for  Eqs.  (3),  (4),  and  (5),  substituting  a  =  k/pCp,  and  multiplying  both 
sides  by  pVCp  yield 


pVCp 

At 


(Tq -T0) 


kAx 

Ax- 


(T,  -  T0) 


+ 


kAx 

Ax+ 


(T2-T0) 


kAy  kAy 

3-T0>  *^<T4-T0> 

+  if!  <t5-To>  *^<t6-t0>  +  Q. 

where  Ax  =  Ay  Az,  Ay  =  Ax  Az,  Az  =  Ax  Ay  and  Q  =  pVCpS. 
x,  y,  and  z  correspond  to  the  coordinates  of  Fig.  2a. 


.'3) 


Fig.  2— Network  of  a  three-dimensional  system 
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The  numbering  system  corresponds  to  a  three-dimensional  network,  shown  in  Fig.  2b. 
It  should  be  obvious  that  the  network  capacitance  value  is  pVCp,  that  the  Gj  value  is 
kAx/Ax--,  etc.  Equation  (6)  may  then  Ik?  written  as 

C0(To  -  T0)/At  =  G,(Tj  T0)  +  G2(T2-T0)  +  G3(T3-T0)  +  G.,(T,  -  T0) 

+  G5(T5  -  T0)  +  Gg(Tg  -  T0)  +  Qy  ,  (7) 

or  in  engineering  terminology,  the  rate  ot  change  of  temperature  with  respect  to  time  is 
proportional  to  the  summation  of  heat  flows  into  the  node. 


It  should  lx?  noted  that  Fig.  2  is  essentially  superimposed  on  a  lumped-parameter  | 

cube  of  a  physical  system  and  is  the  network  representation  of  Eq.  (1).  Since  Eq.  (7)  is  | 

written  in  explicit  form,  only  one  unknown  (Ty)  exists  and  all  of  the  information  neces-  5 

sary  for  its  solution  is  contained  in  the  network  description.  If  it  had  been  formulated 
implicitly,  it  would  have  to  be  solved  in  a  simultaneous  manner.  No  matter  what  method  J 

of  solution  is  requested  of  the  program,  the  information  necessary  has  been  conveyed  by  !; 

the  network  description.  When  an  implicit  set  is  used  with  0  >  0,  the  energy  flows  based  \ 

on  old  temperatures  are  added  to  the  Q  term  and  the  equations  are  then  treated  in  the  | 

same  manner  as  for  0  =  0;  \ 


aV2T  +  S  -  0  .  (8> 

The  solution  of  Pnisson’s  equation  (8)  is  the  solution  utilized  for  steady  state  analy¬ 
sis.  It  is  extremely  important  because  virtually  all  of  the  unconditionally  stable  implicit 
methods  reduce  to  it.  If  Eq.  (7)  had  all  the  right  side  values  primed  and  the  left  side 
was  subtracted  from  both  sides,  we  could  think  of  C0/At  as  a  G0  term  and  T0  (old) 
would  then  become  a  boundary  node.  In  a  manner  of  speaking,  the  capacitor  we  look  at 
in  three  dimensions  becomes  a  conductor  in  four  dimensions.  We  could  draw  a  four- 
dimensional  network,  but  since  there  is  no  feedback  *n  time  it  is  senseless  to  take  more 
than  one  time  step  at  a  time.  However,  various  time-space  transformations  can  be  uti¬ 
lized  such  that  a  one-dimensional  “transient”  yields  the  solution  to  a  two-dimensional 
steady  state  problem,  etc.  This  is  analogous  to  the  “Particle  in  Cell”  method  developed 
in  the  nuclear  field  for  following  shock-wave  propagation. 


Iterative  Techniques 

Now  that  we  have  discussed  the  correlation  between  the  physical,  network,  and 
mathematical  models,  let’s  investigate  the  commonality  of  the  various  methods  of  solu¬ 
tion.  By  describing  the  network  of  Fig.  1  to  the  program,  we  have  supplied  it  with  five 
temperatures,  five  capacitors,  five  sources  (four  not  specified  and  therefore  zero),  four 
conductors,  and  the  adjoining  node  numbers  of  the  conductors.  An  explicit  formulation 
such  as  Eq.  (6)  has  only  one  unknown.  Its  solution  is  easily  obtainable  as  long  as  any 
ass<  dated  stability  criteria  are  continuously  satisfied.  A  more  interesting  formulation 
would  be  a  set  of  implicit  equations  as  follows: 
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(T,'-T1)Ci/At  =  Q\  +  G^-T,') 

(To  ~  T2)C2/At  =  Q'2  +  G^Tj'  -T2)  +  G2(T3  -  T2) 

(T3  -  T3)C3/At  =  Q;}  +  G2(T'  -  T3)  +  G3(T4'  -  T3)  (9) 

(T4-T4)C4/At  =  Q'4  +  G3(T^-T4)  +  G4(T£-T4') 

(T5  -  T5)C5/At  =  Q's  +  G4(T4  -  T5) . 


If  the  above  had  been  formulated  as  partly  explicit  and  implicit,  the  known  explicit  por¬ 
tion  would  have  been  calculated  and  added  to  the  Q  terms,  then  the  0  factor  would  have 
been  divided  into  the  terms  and  multiplied  by  the  At  term. 

If  we  divide  the  At  term  into  the  C  terms  and  indicate  this  by  priming  C,  we  can 
reformulate  Eq.  (9)  as  follows: 


(c;+G,)  T ;  -  Q\  +  c;tj  ♦  g 

(C^  +  G,  +  G2)T2  =  Q'o  +  C2T2  +  G,Tj'  +  G2T3 

(C3  +  G2  +  G3)T3  =  Q3  f  C3T3  -  g2t2  +  G3T4'  (10) 

(C4  +  g3  +  G4  )T4  -  Q'  +  c't4  +  g3t3  +  g4t^ 


<c;  +  g4) 


T5  =  Qs  +  c‘5t5  +  g4t4'  . 


This  equation  can  be  generalized  as 


t;  - 


Cj'Tj  +  xg3t3'  Hr  q\ 
c;  +  eg. 


(id 


where  the  subscript  a  indicates  connection  to  adjoining  nodes.  A  C'  value  of  zero  yields 
the  standard  steady  state  equation,  the  conductor  weighted  mean  of  all  the  surrounding 
nodes.  We  see  here  that  the  C'  can  be  thought  of  as  a  conductor  to  the  old  temperature 
value  and  therefore  Eq.  (11),  although  utilized  to  obtain  transient  solutions,  can  be  con¬ 
sidered  as  a  steady  state  equation  in  four  dimensions.  By  rewriting  Eqs.  (10)  in  the  form 
of  Eq.  (11)  we  are  in  a  position  to  discuss  iterative  techniques.  By  assuming  all  old 
values  on  the  right  hand  side  of  Eq.  (10),  we  couid  calculate  a  new  set  of  temperatures 
on  the  left  which,  although  wrong,  are  closer  to  the  correct  answer.  This  single  set  of 
calculations  is  termed  an  iteration.  By  replacing  all  of  the  old  temperatures  with  those 
just  calculated,  we  can  then  perform  another  iteration.  This  process  is  called  “block” 
iteration.  A  faster  method  is  to  utilize  only  one  location  for  each  temperature.  This  way, 
the  newest  temperature  available  is  always  utilized.  This  method  is  termed  “successive 
point”  iteration  and  is  generally  25%  faster  than  “block”  iteration.  The  iterative  process 
is  continued  a  fixed  (set  by  user)  number  of  times  or  until  the  maximum  absolute  differ¬ 
ence  between  the  r.cw  and  old  temperature  values  is  less  than  some  prespecified  value 
(set  by  user). 
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Although  the  above  operations  are  similar  to  a  relaxation  procedure,  there  is  a  slight 
difference.  We  are  perfoiming  a  set  of  calculations  in  a  fixed  sequence.  A  relaxation 
procedure  would  continuously  seek  the  node  with  the  maximum  temperature  difference 
between  old  and  new  and  calculate  it.  Programmingwise,  as  much  work  is  required  in 
the  seeking  operation,  which  must  be  consecutive,  as  in  the  calculation.  For  this  reason 
it  would  be  wasteful  to  code  a  true  relaxation  method. 

In  addition  to  the  iterative  approach,  several  solution  subroutines  utilize  an  accelera¬ 
tion  feature  and/or  a  different  convergence  criteria.  Once  it  can  be  determined  that  the 
temperatures  are  approaching  the  steady  state  value,  an  extrapolation  is  applied  in  an 
attempt  to  accelerate  convergence.  This  convergence  criterion  is  the  maximum  absolute 
temperature  change  allowed  between  iterations.  This  criteria,  however,  is  generally  one 
sided  and  any  associated  errors  are  accumulative.  In  o.dei  tc  obtain  greater  accuracy, 
some  subroutines  are  coded  to  perform  an  energy  balance  on  the  entire  system  (a  type 
of  Green’s  function)  and  apply  successively  more  severe  convergence  criteria  until  the 
system  energy  balance  (energy  in  minus  energy  out)  is  within  some  prespecified  tolerance. 


Pseudo-Compute  Sequence 

A  set  of  simultaneous  equations  such  as  Eqs.  (10)  is  quite  often  treated  by  matrix 
methods  and  formulated  as  follows: 


where 


{A}  = 


and 


(A)  {T'J  =  {B}  , 


(Cl  +  Gj ) 

"G, 

0 

0 

0 


(Cg  +  Gj  +  G2) 
"Gn 
0 
0 


-Go 


(C3  ^2  +  G3; 


0 

0 

-G3 

(c;  +  g3  +  g4) 

-G4 


0 

0 

0 

-G, 


(12) 


(Gg  +  G4) 


(13) 


t;' 

Q'i  +  CjTj 

Q'2  +  C2T2 

n 

{B}  =  - 

Q3  +  C3T3 

T4' 

Q4  +  C^T4 

Ji 

Qs  +  c.;t5 

The  inverse  of  (A]  is  then  calculated  and  the  solution  obtained  by  matrix  multiplication; 
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{T'}  -  [A]  _1  {B}  .  (14) 

It  should  be  noted  that  the  one-dimensional  problem  has  no  more  than  three  finite 
values  in  any  row  or  column  of  the  coefficient  matrix  [A].  A  three-dimensional  problem 
would  generally  ha'e  no  more  than  seven  finite  values  in  any  row  or  column.  It  is  easy 
to  see  that  a  1000-node,  three-dimensional  problem  would  require  one  million  data  loca¬ 
tions,  of  which  approximately  993,000  would  contain  zero.  The  inverse  might  require 
an  additional  one  million  data  locations.  Aside  from  exceeding  computer  core  area,  the 
computer  time  required  to  calculate  the  inverse  is  proportional  to  the  cube  of  the  prob¬ 
lem  size,  and  large  problems  soon  become  uneconomical  to  solve. 

The  explicit  and  iterative  implicit  methods  previously  discussed  are  well  suited  for 
optimizing  the  data  storage  area  required.  Note  the  adjoining  node  numbers  associated 
with  the  conductors  of  Fig.  1; 

1.1.2  -  G1  between  nodes  1  and  2 

2.2.3  -*•  G2  between  nodes  2  and  3 

(15) 

3.3.4  -*  G3  between  nodes  3  and  4 

4.4.5  -»  G4  between  nodes  4  and  5 . 

Note  also  the  row  and  column  position  of  conductor  values  off  the  main  diagonal  in  the 
(A]  coefficient  matrix,  Eq.  (13).  By  retaining  the  adjoining  node  numbers  for  each  con¬ 
ductor  we  are  able  to  identify  its  element  position  in  the  coefficient  matrix.  As  a  con¬ 
sequence  we  need  store  only  the  finite  values.  The  main  diagonal  term  in  a  composite  of 
the  node  capacitance  and  conductor  values  off  the  main  diagonal. 

The  program  operates  on  the  adjoining  node  numbers  to  form  what  is  termed  the 
pseudo-compute  sequence  (PCS).  The  nodes  are  to  be  calculated  sequentially  in  ascend¬ 
ing  order,  so  the  adjoining  nodes  are  searched  until  the  number  1  is  found.  When  this 
occurs  the  conductor  number  and  the  adjoining  node  number  are  stored  as  a  doublet 
value.  The  search  is  continued  until  all  nodes  one  are  located  and  the  conductor  number 
for  the  last  receives  a  minus  sign.  The  process  is  then  continued  for  node  two,  etc.,  until 
all  the  node  numbers  have  been  processed.  The  pseudo-compute  sequence  formed  (LPCS) 
is  shown  below  left.  A  slight  variation  to  this  operation  is  to  place  a  minus  sign  on  the 
original  other  adjoining  node  number  so  that  it  is  not  recognized  when  it  is  searched  for. 
The  resulting  pseudo-compute  sequence  thus  formed  (SPCS)  is  shown  below  right. 

LPCS  SPCS 

-1,2  -1,2 

1.1  "2,3 

-2,3  -3,4 

2.2  -4,5 

-3,4  -0,0 

3.3 
-4,5 
-4,4 
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The  above  pseudo-compute  sequences  are  termed  long  (LPCS)  and  short  (SPCS),  respec¬ 
tively.  By  starting  at  the  top  of  the  pseudo-compute  sequence,  we  are  operating  on  node 
1.  The  two  values  identify  the  conductor  into  the  node  (the  position  of  the  conductor 
value  in  an  array  of  conductor  values)  and  the  adjoining  node  (the  position  of  the  tem¬ 
perature,  capacitor,  and  source  values  in  arrays  of  temperature,  capacitor,  and  source 
values,  respectively).  The  node  being  operated  on  starts  as  one  and  is  advanced  by  one 
each  time  a  negative  conductor  number  is  passed. 

It  is  easy  to  see  that  the  LPCS  identifies  the  element  position  and  value  locations  of 
all  the  off-diagonal  elements  of  the  row  being  operated  on.  It  takes  complete  advantage 
of  the  sparsity  of  the  coefficient  matrix.  It  is  well  suited  for  “successive  point”  iteration 
of  the  implicit  equations  because  all  elements  in  a  row  are  identified.  When  a  row  is 
processed  and  the  new  T  value  obtained,  the  new  T  can  then  be  used  in  the  calculation 
procedure  of  succeeding  rows. 

The  SPCS  identifies  each  conductor  only  once  and  in  this  manner  takes  advantage 
of  the  symmetry  of  the  coefficient  matrix  as  well  as  the  sparsity.  It  is  well  suited  for 
explicit  methods  of  solution.  The  node  being  operated  on  and  the  adjoining  node  num¬ 
ber  reveal  their  temperature  value  locations  and  their  source  value  locations.  The  explicit 
solution  subroutines  calculate  the  energy  flow  through  the  conductor,  add  it  to  the  source 
location  of  the  node  being  worked  on,  and  subtract  it  from  the  source  location  for  the 
adjoining  node.  However,  if  the  short  pseudo-compute  sequence  were  utilized  for  im¬ 
plicit  methods  of  solution,  they  would  require  the  use  of  slower  “block”  iterative  proce¬ 
dures.  The  succeeding  rows  do  not  have  all  of  the  elements  defined,  and  the  energy  rates 
passed  ahead  were  based  on  old  temperature  values. 


Data  Logistics 

The  pseudo-compute  sequence  formulated  as  shown  above  allow  the  program  to  store 
only  the  finite  values  in  the  coefficient  matrix,  thereby  taking  advantage  of  its  sparsity. 

In  addition,  the  SPCS  takes  advantage  of  any  symmetry  which  may  exist.  Multiply- 
connected  conductors  which  will  be  covered  in  the  next  section  allow  the  user  to  take 
advantage  of  similarity  as  well.  The  foregoing  is  fairly  easy  to  follow,  especially  if  the 
nodes  and  conductors  start  with  the  number  1  and  continue  sequentially  with  no  missing 
numbers.  This  restriction  is  too  limiting  for  general  use  on  large  network  models.  To 
overcome  this  restrict'  jr.  the  program  assigns  relative  numbers  (sequential  and  ascending) 
to  the  incoming  node  data,  conductor  data,  constants  data,  and  array  data  in  the  order 
received.  Any  numbers  missing  in  the  actual  numbering  system  set  up  by  the  user  are 
packed  out,  thereby  requiring  only  as  much  core  space  as  is  actually  necessary. 

All  solution  (Execution)  subroutines  require  three  locations  for  diffusion  node  data 
(temperature,  capacitance,  and  source)  and  one  location  for  each  conductor  value.  They 
also  may  require  from  zero  to  three  extra  locations  per  node  for  intermediate  data 
storage.  Each  node  in  a  three-dimensional  network  has  essentially  six  conductors  con¬ 
nected  to  it  but  only  three  are  unique;  i.e,,  each  additional  node  requires  only  three 
more  conductors.  Hence,  each  node  in  a  three-dimensional  system  requires  from  six  to 
nine  storage  locations  for  data  values  (temperature,  capacitance,  source,  three  conductors, 
and  up  to  three  intermediate  locations).  The  two  integer  values  that  make  up  a  doublet 
of  the  PCS  are  packed  into  a  single  core  location.  Hence,  for  a  three-dimensional  network, 
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each  node  requires  approximately  three  locations  for  data  addressing  for  the  short  and 
six  locations  for  the  long  pseudo-compute  sequence.  The  number  of  core  locations  re¬ 
quired  for  node  can  vary  from  9  to  1 5. 

The  program  requires  the  user  to  allocate  an  array  o?  data  locations  to  be  used  for 
intermediate  data  storage  and  initialize  array  start  and  length  indicators.  Each  subroutine 
that  requires  intermediate  storage  area  has  access  to  this  array  and  the  start  and  length 
indicators.  They  check  to  see  that  there  is  sufficient  space,  update  the  start  and  length 
indicators,  and  continue  with  their  operations.  If  they  call  upon  another  subroutine  re¬ 
quiring  intermediate  storage,  the  secondary  subroutine  repeats  the  check  and  update 
process.  Whenever  any  subroutine  terminates  its  operations  it  returns  tb-  start  and  length 
indicators  to  their  entry  values.  This  process  is  termed  dynamic  storage  allocation  and 
allows  subroutines  to  share  a  common  working  area. 


Order  of  Computation 

A  problem  data  deck  consists  of  data  and  operations  blocks  which  are  preprocessed 
by  CINDA  and  passed  on  to  the  system  FORTRAN  compiler.  The  operations  blocks  are 
named  EXECUTION,  VARIABLES  1,  VARIABLES  2  and  OUTPUT  CALLS.  The 
FORTRAN  compiler  constructs  these  blocks  as  individual  subroutines  with  the  entry 
names  EXECTN,  VARBL1,  VARBL2  and  OUTCAL,  respectively.  After  a  successful  com¬ 
pilation,  control  is  passed  to  the  EXECTN  subroutine.  Therefore,  the  order  of  computa¬ 
tion  depends  on  the  sequence  of  subroutine  calls  placed  in  the  EXECUTION  block  by 
the  program  user.  No  other  operations  blocks  are  performed  unless  called  upon  by  the 
user  either  directly  by  name  or  indirectly  from  some  subroutine  which  int  ynally  calls 
upon  them.  The  network  execution  subroutines  listed  on  p.  internally  call  upon 
VARBL1,  VARBL2,  and  OUTCAL.  Their  internal  order  of  computation  is  quite  simi¬ 
lar,  the  primary  difference  being  the  analytical  method  by  which  they  solve  the  network. 
Figure  3  represents  a  flow  diagram  of  all  the  network  solution  subroutines;  the  subroutine 
writeups  contain  the  comparisons  made  at  the  various  check  points  and  the  routings  taken. 


Systems  Programming 

CINDA  is  actually  an  operating  system  rather  than  an  applications  program.  Two 
programs  are  run  and  executed,  the  second  program  being  the  product  of  the  first.  The 
initial  program,  the  preprocessor ,  operates  in  an  integral  fashion  with  a  large  library  of 
assorted  subroutines  which  can  be  called  in  any  sequence  desired.  It  reads  all  of  the  in¬ 
put  data,  packs  them,  assigns  relative  numbers,  forms  the  pseudo-compute  sequence,  and 
writes  the  data  on  two  different  peripheral  units.  One  unit  contains  FORTRAN  source 
language  generated  from  the  operations  blocks,  with  all  of  the  data  values  dimensioned 
exactly  in  name  COMMON.  This  program,  the  processor,  is  then  compiled  and  executed, 
using  as  input  the  data  from  the  first-mentioned  peripheral  unit.  The  FORTRAN  allo¬ 
cator  has  access  to  the  CINDA  subroutine  library  and  loads  only  those  subroutines  re¬ 
ferred  to  by  the  problem  being  processed. 

Due  to  this  type  of  operation,  CINDA  is  extremely  dependent  on  the  system  soft¬ 
ware  supplied.  However,  once  the  program  has  been  made  operational  on  a  particular 
machine,  the  problem  data  deck  prepared  by  the  user  can  be  considered  as  machine 
independent. 
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III.  DATA  INPUT  REQUIREMENTS 

A  CINDA  problem  data  deck  consists  of  both  data  and  instruction  cards.  The  card¬ 
reading  subroutines  for  CINDA  do  not  utilize  a  fixed  format  type  of  input;  but  rnstead 
use  3  free-form  format.  The  type  of  data  is  oesignated  by  a  mnemonic  code  in  columns 
8,  9,  and  10.  This  is  followed  by  the  data  field  which  consists  of  columns  12-80  or  the 
instruction  field  which  consists  of  columns  12-72.  Although  blanks  are  allowed  before 
or  after  numerical  data  they  may  not  be  contained  within.  The  number  1.234  is  fine, 
but  1.  234  will  cause  the  program  to  abort.  The  program  processes  the  problem  data 
into  FORTRAN  common  data  and  reforms  instructions  into  FORTRAN  source  language 
which  are  then  passed  on  to  the  system  FORTRAN  compiler.  Instruction  cards  which 
contain  an  F  in  column  1  are  passed  on  exactly  as  received.  Any  instruction  card  with 
or  without  an  F  in  column  1  may  contain  a  statement  or  sequence  number  in  columns 
2-5  which  is  passed  on  to  and  used  by  the  FORTRAN  compiler. 

The  most  frequently  used  mnemonic  code  is  three  blanks.  The  data  following  this 
blank  mnemonic  code  may  be  one  or  more  integers,  floating-point  numbers  (with  or  with¬ 
out  the  E  exponent  designation)  or  alphanumeric  words  of  up  to  six  characters  each. 

The  reading  of  a  word  or  number  continues  until  a  comma  is  encountered  and  then  the 
next  word  or  number  is  read.  As  many  numbers  or  words  as  desired  may  be  placed  on 
a  card,  but  they  may  not  be  broken  between  cards.  A  new  card  is  equivalent  to  starting 
with  a  comma  and  therefore  no  continuation  designation  is  required.  AH  blanks  are 
ignored  and  reading  continues  until  the  terminal  column  is  reached  or  a  dollar  sign  en¬ 
countered.  Comments  pertinent  to  a  data  card  may  be  placed  after  a  dollar  sign  and  are 
not  processed  by  the  program.  If  sequential  commas  are  encountered,  floating-point  zero 
values  are  placed  between  them. 

The  next  most  frequently  used  code  is  BCD  (for  binary-coded  decimal)  which  must 
be  followed  by  an  integer  1  through  9  in  column  12.  The  integer  designates  the  number 
of  6-character  words  immediately  following  it.  Blanks  are  retained  and  only  the  desig¬ 
nated  number  of  6- character  words  are  read  from  the  card.  The  mnemonic  code  END  is 
utilized  to  designate  the  end  of  a  block  of  input  to  the  program.  The  code  REM  serves 
the  same  function  as  a  FORTRAN  comment  card;  it  is  not  processed  by  the  program 
but  allows  the  user  to  insert  nondata  for  clarification  purposes.  The  special  codes  CGS, 
CGD,  and  GEN  will  be  discussed  later  in  this  section. 

The  data  deck  prepared  by  a  program  user  consists  of  various  input  blocks  contain¬ 
ing  either  data  or  instructions.  A  fixed  sequence  of  block  input  is  required,  and  each 
block  must  start  with  a  BCD  3  header  card  and  terminate  with  an  END  card  (mnemonic 
codes).  Specific  details  about  these  blocks  follows. 


Title  Block 

The  first  card  of  a  problem  data  deck  is  the  title  block  header  card.  It  conveys  in¬ 
formation  to  the  program  as  to  the  type  of  problem,  which  data  blocks  follow,  and  how 
they  should  be  processed.  The  three  options  presently  available  are 
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BCD  3GBNBRAL 
or  BCD  3THERMAL  SPCS 
or  BCD  3THEHMAL  LPCS 

The  GENERAL  indicates  that  a  nornetwork  problem  follows  and  therefore  no  node  or 
conductor  data  is  present.  The  THERMAL  cards  indicate  that  a  conductor-capacitor 
(CG)  network  description  follows  and  that  either  a  short  (SPCS)  or  long  (LPCS)  pseudo¬ 
compute  sequence  should  be  constructed.  The  title  block  header  card  may  be  followed 
by  as  many  BCD  cards  as  desired.  However,  the  first  20  words  (six  characters  each)  are 
retained  by  the  program  and  used  as  a  page  heading  by  the  user-designated  output  rou¬ 
tines.  The  block  must  be  terminated  by  an  END  card  and  is  then  followed  by  node  data 
for  a  CG  network  problem  or  constants  data  for  a  nonnetwork  problem. 


Node  Data  Block 

As  discussed  in  Section  II,  there  are  three  typos  of  nodes;  diffusion,  arithmetic,  and 
boundary.  Diffusion  nodes  are  those  nodes  with  a  positive  capacitance  and  have  the 
ability  to  store  energy.  Their  future  values  are  calculated  by  a  finite  difference  represen¬ 
tation  of  the  diffusion  partial  differential  equation.  Arithmetic  nodes  are  designated  by 
a  negative  capacitance  value;  they  have  no  physical  capacitance  and  are  unable  to  store 
energy.  Their  future  values  are  calculated  by  a  finite  difference  representation  of  Pois¬ 
son’s  partial  differential  equaiion.  This  is  a  steady  state  calculation  which  always  utilizes 
the  latest  diffusion  node  values  available.  Boundary  nodes  are  designated  by  a  minus 
sign  on  the  node  number;  they  refer  to  the  mathematical  boundary,  not  necessarily  the 
physical  boundary.  Their  values  are  not  changed  by  the  network  solution  subroutines 
but  may  be  modified  as  desired  by  the  user. 

A  diffusion  node  causes  three  core  locations  to  be  utilized,  one  each  for  tempera¬ 
ture,  capacitance,  and  a  source  location.  An  arithmetic  node  receives  core  locations  only 
for  temperature  and  source  and  a  boundary  node  receives  only  a  temperature  location. 
The  program  user  is  required  to  group  his  node  data  into  the  above  three  classes  and  sub¬ 
mit  them  in  that  order.  Node  data  input  with  the  three-blank  mnemonic  code  always 
consists  of  three  values— the  integer  node  number  followed  by  the  floating-point  initial 
temperature  and  capacitance  values.  A  negative  capacitance  value  is  used  to  designate  an 
arithmetic  node  while  a  negative  node  number  designates  a  boundary  node.  Although  the 
capacitance  value  of  a  boundary  node  is  meaningless,  it  must  be  included  so  as  to  main¬ 
tain  the  triplet  format. 

All  nodes  are  renumbered  sequentially  (from  one  on)  in  the  order  received.  The 
user  input  number  is  termed  the  actual  node  number  while  the  program  assigned  number 
is  termed  the  relative  node  number.  This  relative  numbering  system  allows  sequential 
packing  of  the  data  and  does  not  require  a  sequential  numbering  system  on  the  part  of 
the  program  user.  It  is  worth  noting  that  the  pseudo-compute  sequence  is  based  on  the 
relative  numbering  system.  Hence,  the  computational  sequence  of  the  nodes  is  identical 
with  their  input  sequence.  If  a  user  desired  to  reorder  the  computations  in  order  to  aid 
boundary  propagation,  he  needs  merely  to  reorder  his  nodal  input  data. 
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The  mnemonic  codes  CGS,  CGD,  and  GEN  may  be  used.  The  CGS  and  CGD  codes 
are  used  when  one  or  two  materials,  respectively,  with  ternperature-varying  properties  are 
to  be  considered.  For  a  single  material  the  node  number  and  initial  temperature  remain 
the  same  but  instead  of  a  capacitance  value,  the  user  may  input  the  starting  location 
(integer  count)  or  a  doublet  array  of  the  temperature-varying  property  followed  by  the 
actual  (literal)  multiplying  factor  value  to  complete  the  calculation  or  a  constants  loca¬ 
tion  containing  it.  For  a  node  consisting  of  two  materials,  the  node  number  and  initial 
temperature  remain  the  same  but  the  user  would  use  two  array  addresses  and  multiplying 
factors  with  a  CGD  code.  These  codes  would  look  as  follows: 

8 

I 

CGS  N#,TI,AI,F1 
or  CGD  N#,TI ,AI ,FI ,A2,F2 

where  N=  is  the  integer  node  number  and  Ti  is  the  floating-point  initial  temperature. 

The  A  arguments  refer  to  doublet  arrays  of  temperature-varying  Cp  or  p* Cp,  and  the  F 
arguments  may  be  or  refer  to  a  constant  location  containing  the  weight  or  volume,  re¬ 
spectively.  The  CGS  code  causes  the  product  of  the  interpolated  value  times  the  F  factor 
to  be  used  as  the  capacitance  value.  The  CGD  code  uses  the  sum  of  the  separate  inter¬ 
polations  times  the  factor  products  as  the  capacitance  value. 

To  input  a  group  of  sequential  nodes,  the  following  code  is  available: 

8 

■f 

GEN  N#,#N,IN,TI,X,Y,Z,W 

where 

N-  is  the  starting  node  number 

=N  is  the  total  number  of  nodes  desired  (integer) 

IN  is  an  increment  for  the  generated  nodes  (integer) 

Ti  is  the  initial  temperature  for  all  nodes, 

and  the  'apacitance  value  is  calculated  as  the  produce  of  X  times  Y  times  Z  times  W.  If 
this  product  is  negative,  arithmetic  nodes  will  be  generated.  If  Ns?  is  negative,  boundary 
nodes  will  be  generated.  A  sample  node  data  block  could  be  as  follows: 

8  12 
I  I 

BCD  3NUDE  DATA 

1  ,80.,  1  .2,2,80.,  I  .3  STV10  DIFFUSION  NUDES 

CGS  3,80. ,AI ,4.63  SSINOLE  MATERIAL  NODE 

CGD  4,80. ,A! ,2.31 ,A2,4.7o  SDOUULE  MATERIAL  NODE 

GEN  5, 10, 1,30. ,4. 63,1.,!., I.  SGENERATb  10  N0DtS,5-l4 
15,80. ,-!  .,16,80.  ,-l.  STViO  ARITHMETIC  NODES 

-18, -460., 0  SONE  BOUNDARY  NODE 

END 
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The  above  does  not  correspond  to  a  problem;  it  just  represents  data  input.  Note  that  the 
nodes  are  input  in  the  order:  diffusion,  arithmetic  and  boundary.  The  factor  portion  of 
the  CGS  and  CGD  codes  may  be  a  literal  (actual  value)  as  shown  or  reference  a  constant’s 
location  containing  the  value.  Either  one  (not  both)  of  the  array  arguments  on  the  CGD 
code  may  be  a  literal  if  the  property  is  constant.  Both  codes  set  up  linear  interpolation 
calls  which  utilize  the  node  temperature  as  the  independent  variable  and  interpolate  a 
dependent  value  which  is  then  multiplied  by  the  factor  to  obtain  the  capacitance  value. 
The  CG1~*  Coiil  causes  two  interpolations  and  multiplications  to  be  performed  and  sums 
the  products  to  obtain  the  capacitance  value.  These  interpolations  are  performed  each 
iteration  during  the  transient  analysis. 

The  GEN  code  expects  values  in  the  following  order;  starting  node  number,  number 
of  nodes  to  be  generated,  an  increment  for  indexing  the  generated  node  numbers,  the 
initial  temperature  for  all  nodes,  and  four  floating  point  numbers,  the  product  of  which 
is  the  capacitance  value. 


Conductor  Data  Block 

Two  basic  types  of  conductors  may  be  used,  regular  or  radiation,  and  either  may 
utilize  temperature-varying  properties  in  calculating  the  conductance  value.  When  utilizing 
the  blank  mnemonic  code  a  regular  conductor  consists  of  the  integer  conductor  number 
followed  by  two  integer  adjoining  node  numbers  and  the  floating-point  conductance 
value.  If  more  than  one  conductor  has  the  same  constant  value,  they  may  share  the 
same  conductor  number  and  value.  This  is  accomplished  by  placing  two  or  more  pairs  of 
integer  adjoining  node  numbers  between  the  conductor  number  and  value.  The  CGS  and 
CGD  mnemonic  codes  may  also  be  utilized  for  conductors.  They  would  appear  as  follows 

8 

1 

CGS  G#,NA,NB,Al ,F1 
or  CGD  G#,NA,NB, Al ,FI ,A2,F2 

where 

G=  is  the  integer  conductor  number 
NA  is  one  adjoining  node  number 
NB  is  the  other  adjoining  node  number. 

The  A  arguments  refer  to  doublet  arrays  of  temperature-varying  thermal  conductivity 
k(T),  and  the  F  arguments  may  be  or  refer  to  a  constant  location  containing  the  cross 
sectional  area  divided  by  path  length. 

For  CGS  with  FI  positive 

G  =  kl(Tm)*Fl,  Tm  =  (Ta  +  Tb)/2.0.  (16) 

For  CGS  with  FI  negative 

G  =  kl(Ta)*|F3 1.  (17) 
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For  CGD 


kl(Ta)*Fl  k2(Tb)*F2 ' 


(18) 


The  CGS  mnemonic  code  may  be  utilized  for  either  regular  or  radiation  conductors. 
The  data  consist  of  the  integer  conductor  number  and  one  pair  only  of  integer  adjoining 
node  numbers,  and  are  followed  by  an  array  address  and  multiplying  factor.  A  regular 
conductor  would  normally  utilize  the  CGS  code  where  the  addressed  array  would  be 
thermal  conductivity  vs  temperature,  and  the  multiplying  factor  would  consist  of  the 
cross-sectional  area  divided  by  path  length.  A  surface  radiation  conductor  would  utilize 
the  CGS  code  for  a  temperature-varying  array  of  emissivity  with  the  multiplying  factor 
being  the  product  of  surface  area  times  the  Stephan-Boltzmann  constant  (F  =  1.0). 


The  CGD  code  may  be  utilized  for  regular  conductors  passing  through  two  mate¬ 
rials.  In  this  case  two  temperature-varying  property  arrays  and  multiplying  factors  are 
input.  Two  conductance  values  are  calculated  and  one  over  the  summation  of  their  in¬ 
verses  is  returned  as  the  conductor  value.  Either  of  the  array  addresses  may  be  a  literal 
if  one  of  the  properties  is  a  constant.  The  GEN  code  is  also  available  for  conductors  and 
is  input  as  follows: 


8 

I 

GEN  G#,rfG,lG,NA,INA,NB,INB,X,Y,Z,W 


where 


Gs  is  the  starting  conductor  number 

~G  is  the  total  number  of  conductors  desired  (integer) 

IG  is  an  increment  for  the  generated  conductors  (integer) 

NA  and  NB  are  initial  adjoining  node  numbers  (integers) 

INA  and  INB  are  increments  for  the  generated  adjoining  nodes  (integers). 


and  all  generated  conductors  receive  the  same  conductance  value  of  X  times  Y  times  Z 
divided  by  W.  A  negative  G=  will  cause  radiation  conductors  to  be  generated. 


The  GEN  code  may  be  used  to  generate  sequential  conductors,  either  radiation  or 
regular.  The  data  consist  of  the  integer  conductor  number,  an  integer  for  the  number 
of  conductors  to  be  generated,  an  integer  increment  for  indexing  the  generated  conduc¬ 
tors,  the  first  integer  adjoining  node  number,  an  integer  increment  for  indexing  the  first 
adjoining  node  number,  the  second  integer  adjoining  node  number,  an  integer  increment 
for  indexing  the  second  adjoining  node  number,  and  finally  four  floating-point  numbers; 
the  product  of  the  first  three  divided  by  the  fourth  is  the  constant  conductance  value. 
For  example: 


8 

l 

GbN  1,2, I, I* I. 2, 1,2. ,2. ,2. ,2. 

GtN  -3,3,0, 1 » I « 10»0, 1 . , I • « I .* 1 .H+15 


is  equivalent  to 
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12 

i 

I  ,1  ,2,4. ,2, 2, 3, 4. 

-3, I, 10,2, 10,3, 10,1 .b-15. 


An  additional  feature  of  the  program  is  the  one-way  conductor.  This  is  a  conductor 
value  which  enters  into  the  temperature  calculation  of  only  one  of  its  adjoining  nodes 
and  is  indicated  by  placing  a  minus  sign  on  the  unaffected  node.  The  CGS,  CGD,  and 
GEN  codes  may  be  used  for  one-way  conductors.  Physically  this  occurs  in  incompres¬ 
sible  fluid  flow,  and  therefore  the  upstream  node  would  receive  the  minus  sign. 

A  program  idiosyncrasy  which  should  be  mentioned  is  that  while  a  single-valued 
conductor  with  as  many  adjoining  node  pairs  as  desired  may  be  used,  extending  several 
cards  if  necessary,  an  adjoining  node  pair  must  not  be  split  between  cards.  In  addition, 
the  CGS,  CGD,  and  GEN  card  may  have  more  than  one  set  of  data  on  a  card,  but  a  set 
of  data  may  not  be  broken  between  cards.  All  regular  conductors  must  be  entered  prior 
to  any  radiation  conductors.  The  following  is  illustrative  of  the  various  conductor  input 
options. 

8 

1 

BCD 


CGS 

CGD 

GEN 

CGS 

GhN 

END 


3C0NDUCT0R  DATA 

1 , 1 ,2,1 .2, 2, 2,3, I .7 

3, 3, 4, 4, 5, 5,o, 1 *5 

4, -7, d, -8, 9, 7. 6 

5,10, 1 1 ,A3,4.6 

A, 12, 13, A3, 4. 1 ,A4,7.A 

7,3,l,l,l,9,l,l.6,4.0,l.,l. 

-10,1,99, I . E- 1  5 

-II  ,2 ,99, A5 , 1 . E-l 4 

-12,4, 1 ,3, I ,99,0, 1 .E-I4,l . , 1 . , 


STWU  REGULAR  CONDUCTORS 
STRIPLE  PLACED  CONDUCTOR 
SDOUBLE  PLACED  ONE-WAY  COND. 
SVARI ABLE  CONDUCTOR,  SINGLE 
SVARI ABLE  CONDUCTOR,  DOUBLE 
SGENERATE  THREE  CONDUCTORS 
S RAD  I  AT I  ON  CONDJCTOR 
SVARI ABLE  EMISSIVITY  RADIATION 
.  SGENERATE  FOUR  RADIATION  COND. 


The  first  GEN  card  is  equivalent  to  the  following: 

12 

1 

7,1,9,6.4,8,2,10,6.4,9,3, 11, A. 4 

and  the  second  GEN  card  is  equivalent  to 

12 

1 

-12,3,99, 1 .E-l 4, -13, 4, 99, I .E-l 4 
-14,5,99, 1 .E-l 4, -15, 5,99, 1. E-l 4 


If  the  second  GEN  card  had  incremented  the  conductor  number  by  zero,  it  would  have 
been  equivalent  to 
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-12,3,99,4,99,5,99,6,99,1  .E-l  4 

Once  the  node  and  conductor  data  have  been  read  bj  the  program,  construction  of 
the  pseudo-compute  sequence  is  performed.  Any  errors  encountered  cause  an  appropriate 
error  message  to  be  printed  and  a  “do  not  execute”  switch  to  be  set.  However,  the  pro¬ 
gram  will  continue  to  process  input  data  and  attempt  to  discover  any  and  all  recognizable 
errors.  Items  checked  for  are  no  duplicate  node  or  conductor  numbers,  all  conductor  ad¬ 
joining  nodes  must  have  been  specified  in  node  data,  and  all  diffusion  and  arithmetic 
nodes  must  have  at  least  one  conductor  into  them.  A  missing  comma  will  dislocate  the 
data  input  sequence  causing  pages  of  error  messages.  If  over  200  error  messages  are 
printed,  the  program  gives  up  and  immediately  terminates. 


Constants  Data  Block 


Constants  data  are  always  input  as  doublets,  the  constant  name  or  number  followed 
by  its  value.  They  are  divided  into  two  types,  control  constants  and  user  constants,  and 
may  be  intermingled  within  the  block.  Control  constants  (==  50)  have  alphanumeric 
names  while  user  constants  receive  a  number.  User  constants  are  simply  data  storage 
locations  which  may  contain  integers,  floating-point  numbers,  or  up  to  6-character  alpha¬ 
numeric  words.  It  is  up  to  the  program  user  to  place  data  in  user  constant  locations  as 
needed  and  supply  the  location  addresses  to  subroutines  as  arguments. 

Control  constant  values  are  communicated  through  program  COMMON  to  specific 
subroutines  which  require  them.  However,  any  control  constant  name  desired  can  be 
used  as  a  subroutine  argument.  Wherever  possible,  control  constant  values  not  specified 
are  set  to  some  acceptable  value.  If  a  required  control  constant  value  is  not  specified,  an 
appropriate  error  message  is  printed  and  the  program  terminated.  It  is  up  to  the  user  to 
check  the  writeups  of  subroutines  he  is  using  to  determine  control  constant  requirements. 
A  list  of  control  constant  names  and  brief  description  of  each  follows;  check  subroutine 
writeups  for  exact  usage. 

ARLXCA  The  maximum  arithmetic  relaxation  change  allowed. 

ARLXCC  The  maximum  arithmetic  relaxation  change  calculated. 

ATMPCA  The  maximum  arithmetic  temperature  change  allowed. 

ATMPCC  The  maximum  arithmetic  temperature  change  calculated. 

BACKUP  If  nonzero,  the  time  step  just  done  is  erased  and  redone. 

BALENG  User-specified  system  energy  balance  to  be  maintained. 

CSGFAC  Stability  criteria  multiplication/division  factor. 

CSGMAX  Maximum  stability  criteria  for  the  network.  1  .r.yr.  ,  . 

CSGMIN  Minimum  stability  criteria  for  the  network.  J  '  J  maX  an  m> 

CSGRAL  Stability  criteria  range  allowed. 

CSGRCL  Stability  criteria  range  calculated. 

DAMPA  Arithmetic  node  damping  factor. 

DAMPD  Diffusion  node  damping  factor. 

DRLXCA  The  maximum  diffusion  relaxation  change  allowed. 

DRLXCC  The  maximum  diffusion  relaxation  change  calculated. 

DTIMEH  Highest  time  step  allowed  (maximum). 

DTIMEI  Input  time  step  for  implicit  solutions. 


»£*.  v 


MARY  E.  GEALY 


DTIMEL 

DTIMEU 

DTMPCA 

DTMPCC 

ENGBAL 

IDCNT 

LAX FAC 

LINECT 

L00PCT 

NL00P 

0PE1TR 

0UTPUT 

PAGECT 

TIMEM 

TIMEN 

TIMENO 

TIME0 


Lowest  time  step  allowed  (minimum). 

Time  step  used  for  all  transient  network  problems. 

The  maximum  diffusion  temperature  change  allowed. 

The  maximum  diffusion  temperature  change  calculated. 

The  calculated  energy  balance  of  the  system. 

A  counter  for  STOREP  and  RECALL  identification  (integer). 

Number  of  iterations  before  radiation  conductor  is  linear  (integer). 

A  line  counter  location  for  program  output. 

Program  count  of  iteration  loops  performed  (integer). 

User  input  number  of  iteration  loops  desired  (integer). 

Causes  output  each  iteration  if  set  nonzero. 

Time  interval  for  activating  OUTPUT  CALLS. 

A  page  counter  location  for  program  output. 

Mean  time  for  the  computation  interval. 

New  time  at  the  end  of  the  computation  interval. 

Problem  stop  time  for  transient  analysis. 

Old  time  at  the  start  of  the  computation  interval,  also  used  as  problem 
start  time,  may  be  negative. 


1TEST,JT£ST,KTEST,LTEST,MTEST  are  dummy  control  constants  with  integer 
names. 

RTEST,STEST,TTEST,UTEST,VTEST  are  dummy  control  constants  with  noninteger 
names. 


The  following  is  representative  of  a  constants  data  block: 


8 

1 

BCD  3CQNSTANTS  °VTA 

T I MEND, I O.C,UUTPUT, 1 .0  $  CONTROL  CONSTANTS 

1,10,2,3,3,7,4,8  SlNTbGERS 

5 , 1  .  ,6 , 1  .  E3 ,7,  1  .E-3  S  FI.  OAT  I  NO  POINT 

8, TEMP, 9, ALPHA  S ALPHANUMERIC 

END 


Array  Data  Block 

Array  data  arc  exceedingly  simple  to  enter.  The  user  inputs  an  array  number,  se¬ 
quentially  lists  his  information,  and  terminates  it  with  an  END  (data  END,  not  mnemonic). 
Numerous  subroutines  (interpolation,  matrix,  etc.)  require  that  the  exact  number  of  values 
in  an  array  be  specified  as  an  integer.  In  order  to  reduce  the  number  of  subroutine  argu¬ 
ments  and  chance  of  error,  the  CINDA  preprocessor  counts  the  number  of  values  in  an 
array  and  supplies  this  integer  count  as  the  first  value  in  the  array.  The  writeup  of  any 
subroutine  whose  array  arguments  require  the  array  integer  count  will  list  the  array  argu¬ 
ment  as  A(IC).  Subroutines  whose  array  arguments  require  the  first  data  value  rather 
than  the  integer  count  will  list  the  array  argument  as  A(DV).  When  a  user  inputs  the 
array  number  as  positive,  the  integer  count  is  calculated  by  the  preprocessor  and  supplied 
as  the  first  value  in  the  array.  For  example, 
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1 

I  ,  I  .6, 2. 4, 3. 8, END 

The  above  array  (Array  1)  contains  three  data  values  and  was  input  as  a  positive 
array.  By  addressing  A1  as  a  subroutine  argument  the  integer  count  3  would  be  the  first 
value  followed  by  1.6, 2.4  and  3.S.  If  the  user  wanted  the  1.6  data  value  to  be  addressed 
the  argument  should  be  Al+1.  The  user  has  the  option  of  placing  a  minus  sign  on  the 
input  array  number.  In  this  event  the  integer  count  of  data  values  in  the  array  is  not 
calculated  or  stored  and  addressing  the  array  as  A1  obtains  the  first  data  value.  For 
example: 

12 

1 

-2, I .6, 2. 4, 3. 8, END 

Entering  the  argument  A2  would  address  the  1.6  data  value;  the  integer  count  is  not 
available.  The  following  is  an  example  of  various  types  of  arrays. 
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BCD  3ARRAY  DATA 

1 ,1 .6,2.4, 3. 8, END  SFLOATING  POINT  NUMBERS 

2, TEMP  I ,TEMP2,END  S ALPHANUMERIC 

3  S ALPHANUMERIC 

BCD  3TEMPERATURE  STUDY 

END 

-4, SPACE, I  00, END  SSPACE  OPTION 

END 

Two  types  of  alphanumeric  inputs  are  shewn  above.  The  first  allows  each  word  to 
be  separated  by  a  comma,  requires  each  word  to  start  with  a  letter,  and  does  not  allow 
the  use  of  blanks.  The  second  requires  use  of  the  BCD  mnemonic  code  and  the  integer 
word  count.  It  allows  use  of  letters,  numbers,  or  characters  anywhere  and  retains  blanks. 
The  space  option  is  an  easy  way  for  the  user  to  specify  a  large  number  of  locations 
which  are  initialized  by  the  preprocessor  ~s  floating-point  zeros.  The  space  option  re¬ 
quires  the  word  SPACE  followed  by  the  lumber  of  locations  to  be  initialized  It  may 
be  used  anywhere  in  an  array  and  as  ma’iy  times  as  desired  as  long  as  total  available  core 
space  is  not  exceeded. 


Program  Control 

Aside  from  the  title  block,  there  are  either  two  or  four  data  blocks  depending  upon 
whether  the  problem  is  GENERAL  cr<  THERMAL,  respectively.  No  matter  which,  there 
are  also  four  operations  blocks  entitled  EXECUTION,  VARIABLES  1,  VARIABLES  2, 
and  OUTPUT  CALLS.  The  operations  or  instructions  called  for  in  these  blocks  determine 
the  program  control.  They  are  preprocessed  by  CINDA  and  passed  on  to  the  system 
FORTRAN  compiler  as  four  separate  subroutines  entitled  EXECTN,  VARBL1,  VARBL2, 
and  OUTCAL,  respectively.  When  the  FORTRAN  compilation  is  successfully  completed, 
control  is  passed  to  the  EXECTN  subroutine  which  sequentially  performs  the  'perations 
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in  the  same  order  as  entered  by  the  user  in  the  EXECUTION  block.  None  of  the  opera¬ 
tions  specified  in  the  other  three  blocks  will  be  performed  unless  they  are  called  for, 
either  directly  by  name  in  the  EXECUTION  block  or  internally  by  some  other  called-for 
subroutine. 


No  operations  will  be  performed  unless  requested  by  the  user,  and,  no  control  con¬ 
stants  will  be  utilized  unless  some  subroutine  calls  for  them.  Network  solution  subrou¬ 
tines  internally  call  upon  VARBL1,  VARBL2,  and  OUTCAL  (see  Fig.  3).  They  also  use 
numerous  control  constants,  but  their  individual  writeups  in  Section  VIII  must  be  con¬ 
sulted  to  determine  which  ones  and  their  exact  usage.  Network  solution  subroutines  re¬ 
quire  no  arguments  but  most  others  do.  These  arguments  may  be  addresses  which  refer 
to  the  location  of  data  or  they  may  be  literals;  i.e.,  the  actual  data  value.  All  of  the 
inpu',  data  can  be  addressed  by  using  alphanumeric  arguments  of  the  following  form. 


TN  for  the  temperature  location  of  node  N 
CN  for  the  capacitance  location  of  node  N 
QN  for  the  source  location  of  node  N 
GN  for  the  conductance  location  of  conductor  N 
KN  for  the  value  location  of  constant  N 
AN  for  the  starting  location  of  array  N 


*  v 


and  control  constants  utilize  their  individual  names. 


mation 
cards 

3EXECUTI0N  card.  For  example, 


7 

I  i 

DIMENSION  X< 
BCD  3EXECUT i OK 
NDIM  =  100 
NTH  =  0 


t>‘.  '» ...  . 


In  the  DIMENSION  card 
must  be  in  an  15  format 
would  cause  a  v'orking  array 
are  needed,  the  integer  100 
NDIM).  If  no  working  locations 
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F  in  co,umn  1  indicates  to  the  preprocessor  that  the  card  is  FORTRAN  and 
shoulu  be  passed  on  as  received.  This  F  option  allows  the  user  to  program  FORTRAN 
operations  directly  into  the  operations  blocks.  However,  the  C1NDA  arguments  listed 
aoove  are  not  FORTRAN  compatible  with  the  exception  of  the  control  constant  names, 
therefore,  it  is  recommended  that  the  program  user  utilize  CINDA  subroutine  calls 
wherever  possible.  This  is  impossible  however  when  logical  operations  are  required  In 
this  case  it  is  recommended  that  the  user  place  CINDA  data  values  as  needed  into  the 
mailable  dummy  control  constant  names  allowed  for  that  purpose.  Then,  FORTRAN 

C3n  1)6  utUized  with  the  dummy  control  constant  names  as  arguments. 
rUK I KAN  statement  numbers  for  routing  purposes  may  be  placed  in  columns  2-5  on 
any  operations  cards. 

i  o  Jh!,data  fiGld  for  node’  conductor>  constant,  and  array  data  consists  of  columns 
,  80,  However,  the  data  field  of  operations  cards  ends  with  column  12.  In  a  manner 
ot  speakmg,  a  CINDA  subroutine  call  is  a  special  array  and  should  terminate  with  a  data 
h-ND.  In  order  to  simplify  input  for  the  user,  the  operations  read  subroutines  recognize 
two  special  characters;  the  left  and  right  parentheses.  The  left  parenthesis  is  accepted  as 
a  comma,  while  the  right  parenthesis  is  accepted  as  a  comma  followed  by  a  data  END 
I  his  allows  what  would  have  been 


r  • 


ADD.KI ,K2#K3,ENn 


B  .a 

k  t0  f,e  more  aesthetically  formatted  as 

.... 


|  -  'f  S  "  j, 

f. 

h  -  f  •  *  ?•- 

. 

feT'  ’r  ‘ 

tit's- 

fep  . 
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ADD(KI  ,K2,K3) 


which  is  almost  identical  to  a  FORTRAN  subroutine  call. 

Execution  Operations  Block 

An  Execution  operation  block  might  be  as  follows: 


7  12  21  2b 

11  11 

DIMENSION  X(  2b) 

BCD  3 EXECUTION 
NDIM  =  25 
NTH  =  0 

10  TIMEND  =  TIMEND  +  1 .0 
CNFRWD 

STFShP(T20,TTEST) 

IF ( TTEST  .LE.  100.)  GO  TO 
END 


$  EXPLICIT  FORWARD  DIFFERENCING 
$  PLACE  TIO  INTO  DUMMY  CC 
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The  above  indicates  a  transient  thermal  problem  in  which  the  user  desires  to  terminate  the 
analysis  when  the  temperature  at  node  20  exceeds  100°F.  The  problem  must  have  been 
fairly  small  because  only  25  working  locations  were  dimensioned  and  CNFRWD  requires 
one  per  node.  It  does  demonstrate  the  use  of  both  CINDA  calls  and  FORTRAN  opera¬ 
tions  and  that  control  constants  are  referred  to  by  name  in  either.  Another  example 
might  be 

1  8  12  21  25 

i  11  11 

F  DIMENSION  X(  500) 

BCD  3EXECUTI0N 
F  NDIM  =500 

F  NTH  =0 

Cl HDSL  $  STEADY  STATE  (USES  LPCS) 

F  TIMEND  =  10.0 

CNFRWD  STRANSI ENT  ANALYSIS  (USES  SPCS) 

END 

In  this  case  the  user  desires  to  have  a  steady  state  analysis  performed  on  the  network  and 
then  a  transient  analysis  performed  utilizing  the  steady  state  answer  as  initial  conditions. 
However,  the  two-network  solution  subroutines  referred  to  are  incompatible  in  their  PCS 
requirements  and  the  program  would  be  terminated  with  an  appropriate  error  message. 

A  further  example  might  be 

8  12 
1  1 

BCD  3EXECUTI0N 

INVERSES  A I  ,A  2) 

MULT(A2,A3) 

LIST(A2,KI ,17) 

LIST(A3,K2,  17) 

END 

The  above  problem  consists  entirely  of  matrix  operations  and  therefore  is  run  as  a  GEN¬ 
ERAL.  The  subroutines  do  not  require  any  working  space  so  none  have  been  dimensioned. 
Furthermore,  no  reference,  direct  or  indirect,  is  made  to  VARBL1,  VARBL2,  or  OUTCAL, 
and  those  operations  blocks  should  be  empty.  Even  though  they  may  be  empty  or  not 
referred  to,  their  blockheader  and  mnemonic  END  cards  must  still  be  entered. 

There  is  no  end  to  the  variety  of  examples  that  could  be  generated.  In  reality,  the 
program  user  is  actually  programming,  although  it  is  comewhat  disguised  as  data  input. 
However,  the  program  does  simplify  the  task  of  data  logistics,  and  it  automates  data  input 
and  output,  construction  of  the  PCS,  loading  the  subroutine  library,  and  other  systems 
features,  thereby  greatly  lessening  the  programming  knowledge  which  might  otherwise  be 
required  of  a  user. 

A  point  well  worth  considering  is  proper  initialization.  All  instructions  contained  in 
the  other  three  operations  blocks  are  performed  each  iteration  or  on  the  output  interval. 

If  an  operation  being  performed  in  Variables  1  is  utilizing  and  producing  nonchanging 
constants,  it  should  be  placed  in  the  Execution  block  (prior  to  the  network  solution  call)  so 
that  it  will  be  performed  only  once.  Input  arrays  requiring  postinterpolation  multiplication 


$ SEE  MATRIX  SUBROUTINE 
SrtRITEUPS  FOR  OPERATIONS 
S PERFORMED 
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for  units  conversion  only  could  be  prescaled,  thereby  deleting  and  multiplication  process. 
Complex  functions  of  a  single  independent  variable  requiring  several  interpolation  values 
which  are  then  combined  in  a  multiplicative  fashion  can  be  precalculated  vs  the  independ¬ 
ent  variable.  Such  a  precalculated  complex  function  reduces  the  amount  of  work  per¬ 
formed  during  the  transient  analysis.  A  great  many  operations  of  this  type  can  be  per¬ 
formed  in  the  Execution  block  prior  to  call  for  a  transient  analysis.  Also,  output  operations 
to  be  performed  once  the  transient  analysis  is  completed  may  be  placed  directly  into  the 
Execution  block  following  the  transient  network  solution  call. 


Variables  1  Operations  Block 

The  statement  that  this  program  solves  nonlinear  partial  differential  equations  of  the 
diffusion  type  is  not  quite  accurate.  In  reality  the  program  only  solves  linear  equations. 
However,  nonlinearities  are  evaluated  at  each  computation  interval  and  in  this  manner 
generally  yield  acceptable  answers  to  nonlinear  problems.  This  method  is  rrore  properly 
termed  quasilinearization.  The  Variables  1  operation  block  allows  a  point  in  the  computa¬ 
tional  sequence  at  v/hich  the  user  can  specify  the  evaluation  of  nonlinear  network  ele¬ 
ments,  coefficients,  and  boundary  values  (see  Fig.  3).  The  CGS  and  CGD  mnemonic 
codes  utilized  for  node  and  conductor  data  cause  the  construction  of  various  subroutine 
calls  which  are  placed  in  this  block  by  the  CINDA  preprocessor.  The  user  must  specify 
any  additional  subroutine  calls  necessary  to  completely  define  the  network  prior  to  enter¬ 
ing  the  network  solution  phase. 

Prior  to  inclusion  of  the  CGS  and  CGD  mnemonic  codes,  the  Variables  1  operations 
block  primarily  consisted  of  linear  interpolation  subroutine  calls  input  by  the  user  for  the 
evaluation  of  temperature  varying  properties  While  these  linear  interpolation  calls  are 
automated  through  use  of  the  CGS  and  CGD  codes,  it  is  up  to  the  program  user  to 
specify  any  required  bivariate  or  trivariate  interpolations  or  other  functional  evaluations 
necessary.  Just  prior  to  performing  the  Variables  1  operations,  all  network  solution  sub¬ 
routines  zero  out  all  source  locations.  Therefore,  the  user  is  required  to  specify  constant 
as  well  as  variable  or  nonlinear  impressed  sources  in  this  block.  A  Variables  1  operations 
block  could  be  as  follows: 

1  8  12 

i  ;  t 

BCD  3VARI ABLfc'S  I 

STFSEPU0.0/JI7)  SCONSTANT  IMPRESSED  SOURCE 

D1DEG1 (TIMEMtA8,Q!9)  STIME  VARYING  SOURCE 

D2DIWM(T1  8fTIMENt  Al  9, 7.<S3fGI8)  SBIVARIATE  FUNCriON 
F  TTEST=  11.6 

F  IF  (TIMEN.GT.  10.>TTcST  =0.0 

STFSEP(TTEST,027)  ^VARIABLE  SOURCE 

END 

The  first  call  above  places  a  constant  heating  rate  of  10.0  into  the  source  location  of 
node  17.  The  second  call  causes  a  linear  interpolation  to  be  performed  on  array  8  using 
mean  time  as  the  independent  variable  to  obtain  a  time-varying  heating  rate  for  node  19. 

The  third  call  uses  mean  time  and  the  temperature  at  node  18  as  independent  variables 
to  perform  a  bivariate  interpolation.  The  interpolated  answer  is  then  multiplied  by  7.63 


25 


MARY  E.  GEALY 


and  placed  as  the  conductance  value  of  conductor  18.  The  next  two  cards  are  FORTRAN 
and  allow  a  value  of  11.6  to  be  placed  into  control  constant  TTEST  until  TIMEN  exceeds 
10.0,  after  which  a  value  of  0.0  if  placed  into  TTEST.  This  amounts  to  a  single  step  in 
a  “staircase”  function.  The  last  card  places  the  value  from  TTEST  into  the  source  loca¬ 
tion  for  node  27.  Another  sample  Variables  1  block  might  look  as  follows: 

8  12 

4  4 

BCD  3VARI A8LES  1 

BLDARYf A! 2+ 1 ,T I ,T7, T3 ,T4)  $CONSTRUCT  VECTOR 

D 1  DEG  1 (T7, A 19, A 1 3+2)  $  INTERPOLATION 

IRRADE(A7,AI3,AI0,A12)  SIR  RADI OSITY  EXPLICIT 

BRKARY(AI2+I ,Q1 ,Q7,Q3,Q4)  SDISTRIBUTE  Q  RATES 

D! DIWMCTI MEM, A9, 0.35,  TTEST) 

ADD(TTEST,QI ,01 )  S ADD  TWO  RATES 

END 

The  first  call  causes  the  construction  of  an  array  of  four  temperature  values  necessary  as 
input  to  an  infrared  radiosity  subroutine  (third  card).  The  second  call  causes  the  linear 
interpolation  of  a  temperature-varying  property  from  array  19  to  be  placed  into  array 
13  +  2  which  is  the  second  array  argument  for  the  radiosity  call.  This  second  argument 
must  be  an  array  of  surface  emissivities  for  the  surfaces  under  consideration;  therefore 
array  19  must  be  an  array  of  temperature-varying  emissivity.  The  BRKARY  call  takes 
data  values  from  array  12  +  1,  2,  3,  and  4  and  places  them  into  the  source  locations  for 
nodes  1,  7,  3,  and  4,  respectively.  The  fifth  call  performs  linear  interpolation  on  array  9 
using  T1MEM  as  the  independent  variable,  multiplies  the  result  by  0.35  and  places  it  in 
control  constant  TTEST.  This  might  be  a  time-varying  solar  heating  rate  where  0.35  is  the 
solar  absorptivity.  The  ADD  call  adds  this  rate  to  what  is  already  contained  in  the  source 
location  for  node  1.  Each  node  has  one  and  only  one  source  location.  If  a  user  desires 
to  impress  more  than  one  heating  rate  on  a  node,  he  must  sum  the  rates  and  supply  the 
value  to  the  single  source  location  available  per  node. 

The  Variables  1  operations  block  is  the  logical  point  in  the  network  computational 
sequence  for  the  calculation  of  impressed  sources  whether  they  are  due  to  internal  dissi¬ 
pation  of  power  components,  radiation  depositation,  aerodynamic  heating,  or  orbital 
heating.  If  a  desired  subroutine  is  not  available,  the  user  may  always  add  his  own;  data 
communication  is  obtained  through  subroutine  arguments  as  in  any  other  subroutine. 


Variables  2  Operations  Block 

With  regard  to  the  network  solution,  the  Variables  1  operations  may  be  thought  of 
as  presolution  operations  and  the  Variables  2  operations  as  postsolution  operations.  In 
Variables  1  the  network  was  completely  defined  with  respect  to  nonlinear  elements  and 
boundary  conditions.  Variables  2  allows  the  user  to  look  at  the  network  just  solved.  He 
may  meter  and  integrate  flow  rates,  make  corrections  in  order  to  account  for  material 
phase  changes,  or  compare  answers  just  calculated  with  test  data  in  order  to  derive  em¬ 
pirical  relationships.  A  simple  Variables  2  operations  block  might  be  as  follows: 
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8  12 
I  I 

BCD  3VARI ABLES  2 

QMETER1T1  ,T2,GI  ,KI ) 
QINTEG(K2,DTIMEU,K2) 
RDTNQS(T5,T1 ,G8,K3) 

Q I NTEG ( K3 , DTI MEU , K4  > 
ADD(K2,K4,K5) 

END 


SMETER  HEAT  FLOW 
$  INTEGRATE  HEAT  FLOW 
SMETER  RADIATION  FLOW 
S INTEGRATE  RADIANT  FI  3W 


The  first  call  measures  the  heat  flow  from  node  1  to  node  2  through  regular  con¬ 
ductor  1  and  stores  the  result  in  constant  location  1.  The  second  call  performs  a  simple 
integration  with  respect  to  time  and  sums  the  result  into  constants  location  2.  The  third 
call  measures  heat  flow  through  a  radiation  conductor  which  is  then  integrated  by  the 
fourth  call.  The  sum  of  the  two  integrations  is  obtained  by  the  fifth  call.  Another  Vari¬ 
ables  2  operations  block  might  be  as  follows: 

8  12 
I  i 

BCD  3VARI ABLES  2 

ABLATS(AI ,1 .76,K8,A7,T15,C15)  SABLATIVE  ON  N0DE15 

END 


Phase  change  subroutines  such  as  the  above  are  unique  in  that  they  perform  auto¬ 
matic  corrector  operations.  Node  15  has  been  solved  by  the  network  solution  subroutine 
as  though  no  ablative  existed.  The  ABLATS  subroutine  then  corrects  the  temperature  at 
node  15  to  account  for  the  ablative  material.  It  does  this  by  calculating  the  average 
heating  rate  to  node  15  over  the  time  step  just  performed  and  utilizes  it  as  an  inner- 
surface  boundary  condition  for  the  internally  constructed  one-dimensional  network  repre¬ 
sentation  of  the  ablative  material.  The  correctness  of  this  analytical  approach  can  be 
rigorously  substantiated  for  use  with  explicit  network  solution  subroutines.  However, 
when  used  with  large  time  step  implicit  methods  it  yields  a  controlled  instability  and  the 
results  may  be  questionable.  It  is  up  to  the  user  to  determine  the  solution  accuracy  by 
whatever  means  available.  A  more  complicated  Variables  2  operations  block  could  be 
as  follows: 

1  5  8  12 

l  i  I  1 

BCD  3VARI ABLES  2 

D1DEGI (TI MEN* Al 0,K8) 

SUB(T8,K8, TTEST) 

F  I F ( TTEST . LE . 2 . 0 ) GU  TO  10 

MLTPLY(G7,0.99,G7) 

5  STFSEP(-1 .0, BACKUP) 

F  GO  TO  20 

F  10  IF(TTEST.GE.-2.0)  GO  TO  lb 
MLTPLYIG7, I .01,07) 

F  GU  TO  5 

15  0METEH(T8,T15,K9) 

QINTEG(K9,DTIMEU, K1 0) 

F  20  CONTINUE 
END 


SGET  TEST  TcMPERATURE 
SOBTAIM  TEMP  DIFFERENCE 

SREDUCE  CONDUCTANCE 
SSET  BACKUP  NON-ZERO 


S INCREASE  CONDUCTANCE 
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Fig.  •!— 1 Three-dimensional 
network 


This  corresponds  to  a  portion  of  a  network  such  as  shown  in  Fig.  4. 

Array  10  is  a  time-temperature  test  history  of  node  8,  and  node  15  is  a  known 
boundary  reference  temperature.  The  problem  is  to  calculate  the  value  of  conductor  7 
which  will  yield  a  calculated  temperature  at  node  8  that  is  within  ±2.0  degrees  of  the  test 
history.  The  above  Variables  2  operations  will  attempt  to  modify  conductor  7  so  that  it 
will  meet  the  constraints  on  temperature  8.  It  is  quite  “brute  force”  and  unsophisticated. 
However,  the  corrector  operations  are  at  the  discretion  of  the  user.  If  the  tolerances 
were  too  severe  or  the  correction  operations  too  strong,  the  correction  for  one  tolerance 
could  lead  to  dissatisfaction  of  the  other  and  an  impasse  result.  If  the  reference  tempera¬ 
ture  at  node  15  were  incorrect,  possibly  no  value  of  conductor  7  would  satisfy  the  con¬ 
straints.  The  end  result  of  such  a  study  would  be  to  produce  a  plot  of  conductance  7  vs 
time  which  could  be  used  to  derive  an  empirical  relationship  with  other  parameters.  Too 
wide  a  tolerance  would  cause  the  plot  to  resemble  a  staircase  function.  Please  note  that 
either  condition  being  unsatisfied  causes  control  constant  BACKUP  to  be  nonzero  and 
the  iteration  to  be  redone  with  the  corrected  conductor  7  value.  Only  when  all  criteria 
are  met  are  the  metering  and  integration  operations  performed. 


Output  Calls  Operations  Block 

This  operations  block  could  have  been  entitled  Variables  3  but  Output  Calls  seemed 
more  appropriate.  In  it  a  user  may  call  upon  any  desired  subroutine.  However,  its  con¬ 
tents  are  performed  on  the  output  interval  (see  Fig.  3),  so  it  is  only  logical  that  it  would 
primarily  contain  instructions  for  outputting  information.  There  is  a  variety  of  output 
subroutines  offering  the  user  several  format  options.  A  very  simple  Output  Calls  block 
would  be  as  follows: 

8  12 
I  i 

BCD  30UTPUT  CALLS 
PRNTMP 

bND 

The  above  call  will  output  certain  time  control  information  and  the  temperature  of  every 
node  in  the  network  under  consideration.  The  node  temperatures  will  correspond  to  the 
relative  node  numbers  set  up  by  the  preprocessor,  not  the  actual  node  numbers  set  by 
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the  user.  The  preprocessor  lists  out  all  of  the  input  data.  Immediately  after  the 
input  node  data  a  dictionary  of  relative  node  numbers  vs  actual  node  numbers  is 
listed.  By  utilizing  it  a  user  can  correlate  the  relative  node  temperatures  with  his 
actual  numbers. 

The  Output  Calls  will  be  performed  at  problem  start  time  and  on  the  output  inter¬ 
val  until  problem  stop  time  is  reached.  For  example,  a  100-min  transient  analysis 
with  an  output  interval  of  5  min  would  cause  the  Output  Calls  operations  to  be  per¬ 
formed  21  times. 

The  above  data  and  operations  blocks  constitute  a  problem  data  deck  which  must  be 
terminated  by  the  following  card: 

8  12 

i  I 

BCD  3t£ND  UF  DATA 


Parameter  Runs 

Parametric  analyses  which  do  not  involve  network  of  operations  changes  to  the 
original  problem  may  be  performed  on  the  same  computer  run.  Only  data  values  such 
as  output  page  heading,  temperatures,  capacitances,  conductances,  constants,  and  arrays 
may  be  changed.  The  data  change  blocks  must  all  be  specified  whether  changes  occur 
in  the  block  or  not,  and  the  data  input  is  identical  to  the  preceding  discus: !  jn  with  the 
exception  of  conductors.  When  specifying  new  conductances,  the  adjoining  node  infor¬ 
mation  is  deleted;  only  the  conductor  number  and  value  are  required.  The  only  mnemon¬ 
ics  allowed  are  the  three  blanks  and  BCD.  When  changing  an  array,  the  entire  new  array 
must  be  entered  and  be  exactly  the  length  of  its  original.  No  new  arrays  or  numbered 
constants  may  be  defined. 

Two  parametric  run  options  are  available,  INITIAL  and/or  FINAL,  and  they  may 
be  used  several  times  within  the  problem  data  deck.  The  problem  data  deck  as  initially 
entered  is  referred  to  as  the  original  problem.  Any  and  all  INITIAL  parameter  runs  refer 
to  it  exactly  as  it  was  put  in.  The  FINAL  parameter  run  refers  to  the  problem  just  com¬ 
pleted  exactly  as  terminated.  When  two  INITIAL  parameter  runs  are  attached  to  the  end 
of  a  problem  data  deck,  they  both  refer  to  the  original  problem  at  start  time.  However, 
when  two  FINAL  parameter  runs  are  attached  to  the  end  of  a  problem  data  deck,  the 
first  refers  to  the  original  as  terminated,  and  the  second  refers  to  the  first  FINAL  param¬ 
eter  run  as  completed.  The  CINDA  control  cards  necessary  to  specify  a  parameter  run 
are  as  follows: 
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8  12 
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BCD  3 INITIAL  PARAMETERS 
or  BCD  3FINAL  PARAMETERS 
END 

BCD  3N0DE  DATA 
END 

BCD  3 CONDUCTOR  DATA 
END 

BCD  3 CONSTANTS  DATA 
END 

BCD  3ARRAY  DATA 
END 


The  parameter  run  decks  are  inserted  in  the  problem  data  deck  immediately  preceding  the 
BCD  3END  OF  DATA  card.  After  the  BCD  parameter  card,  the  user  may  insert  addi¬ 
tional  BCD  data  to  replace  the  original  problem  output  page  heading.  Parameter  runs 
conserve  machine  time  mainly  because  the  PCS  does  not  have  to  be  reformed.  If  a  user 
desires,  he  may  accomplish  FINAL  parameter  runs  by  calling  the  network  execution  sub¬ 
routine  twice  in  the  Execution  block  and  inserting  the  necessary  calls  to  modify  data  val¬ 
ues  between  them. 


Store  and  Recall  Problem  Options 

The  purpose  of  the  store  and  recall  options  is  to  provide  the  user  with  the  means  to 
interrupt  his  program  at  any  point,  store  the  current  data  values  on  tape,  and  continue 
processing.  While  the  parameter  run  capability  is  useful  for  performing  parametric  analy¬ 
ses  in  the  same  ran,  the  store  and  recall  capability  allows  an  indefinite  time  lapse  between 
parametric  analyses.  In  addition,  long  duration  problems  may  be  broken  into  several 
short  duration  runs.  If  a  parametric  analysis  is  such  that  the  first  portion  of  the  runs  are 
identical,  then  the  problem  can  be  run  for  the  constant  portion,  stored  and  then  recalled 
as  many  times  as  necessary. 

The  store  problem  feature  is  achieved  by  a  user  initiated  subroutine  call  which  is 
as  follows: 

12 

1 

STOREPtKX) 

where  KX  refers  to  a  constant  location  containing  an  alphanumeric  identification  name 
for  the  stored  problem.  The  call  may  be  used  as  many  times  as  desired,  but  each  activa¬ 
tion  must  reference  a  unique  name.  To  allow  the  STOREP  call  to  be  placed  in  a  loop, 
the  programmer  must  use  the  control  constant  IDCNT.  By  incrementing  this  constant 
within  the  loop  (not  to  exceed  999),  a  unique  combination  of  KX  and  IDCNT  will  be 
available  to  identify  the  problem. 

The  recall  problem  feature  is  a  CINDA  preprocessor  option,  which  is  activated  by 
the  following  card: 
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RECALL  AAAAAA  NNN 

where  AAAAAA  is  the  alphanumeric  identification  name  of  the  stored  problem,  and  NNN 
is  the  integer  IDCNT.  (If  IDCNT  was  not  used  with  STOREP,  NNN  =  0).  This  single 
card  replaces  the  blank  card  preceding  the  problem  data  deck  and  must  be  followed  by  a 
BCD  3INITIAL  PARAMETERS  data  deck.  The  stored  problem  identified  will  be  searched 
for  and  brought  into  core  from  the  two  storage  tapes.  Any  data  changes  specified  will  be 
performed  and  the  control  is  passed  to  the  first  subroutine  call  in  the  EXECUTION  block. 

The  problem  is  stored  on  logical  unit  22  and  recalled  from  unit  21;  the  processor 
(tape  40)  must  be  saved  in  a  store  run  and  remounted  on  a  recall  run.  The  user  must 
remember  that  the  recalled  problem  contains  the  STOREP  call.  Because  of  this  feature, 
the  user  has  the  option  whether  or  not  to  store  the  problem  again.  Logical  unit  22  is 
equipped  depending  on  this  option.  Section  VII  should  be  consulted  for  details  concern¬ 
ing  deck  setups,  tape  usage,  and  Job  Request  forms. 


Data  Printing  Option 

At  times,  the  user  may  wish  to  see  what  is  being  stored  on  the  data  tapes  during 
preprocessing  (i.e.,  LUT1,  LUT2,  LUT3,  LUT4,  and  LB3D).  A  printout  will  occur  if  an 
asterisk  is  put  in  column  80  of  the  first  BCD  card  (the  GENERAL,  THERMAL,  INITIAL 
PARAMETERS,  or  FINAL  PARAMETERS  card).  After  each  block  of  the  problem  deck 
is  printed,  there  will  be  a  listing  of  tiie  appropriate  data  (most  of  which  will  be  binary). 


IV.  ERROR  MESSAGES 

Due  to  the  variety  of  subroutines  available  and  the  variable  number  of  arguments 
which  some  of  them  have,  no  check  is  made  to  determine  if  a  subroutine  call  has  the  cor¬ 
rect  number  of  arguments.  An  incorrect  number  of  arguments  on  a  subroutine  call  will 
generally  cause  job  termination  immediately  after  successful  compilation,  usually  without 
any  error  message.  If  the  above  occurs,  the  user  should  closely  check  the  number  of  ar¬ 
guments  for  his  subroutine  calls. 

Numerous  error  messages  can  be  put  out  by  the  preprocessor.  These  error  messages 
are  listed  below  and  are  grouped  according  to  various  preprocessor  functions.  All  error 
messages  are  preceded  by  three  asterisks,  which  have  been  deleted  below.  Self-explanatory 
messages  are  not  enlarged  upon. 


Processing  Data  Blocks— 

DATA  BL0CKS  IN  IMPR0PER  0RDER  0R  ILLEGAL  BL0CK 
DESIGNATION  ENCOUNTERED. 

AN  IMBEDDED  BLANK  HAS  BEEN  ENCOUNTERED  IN  THE  LAST  LINE. 
BLANK  COUNT  0F  10  HAS  BEEN  EXCEEDED. 
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INTEGER  FIELD  EXCEEDS  10. 

REAL  NUMBER  FIELD  EXCEEDS  20. 

ALPHAMERIC  FIELD  EXCEEDS  6. 

MULTIPLE  DECIMAL  P0INTS  HAVE  BEEN  ENC0UNTERED. 

N0DES  MUST  BE  0RDERED  -  DIFFUSI0N,  ARITHMETIC,  B0UNDARY. 

C0NDUCT0RS  MUST  BE  0RDERED  -  REGULAR,  RADIATI0N. 

N0DE  NUMBER,  XXXXX,  IS  THE  DUPLICATE  0F  THE  XXXXXTH  N0DE. 

C0NDUCT0R  NUMBER,  XXXXX,  IS  THE  DUPLICATE  0F  THE  XXXXXTH 
C0NDUCT0R. 

C0NSTANT  NUMBER,  XXXXX,  IS  THE  DUPLICATE  0F  THE  XXXXXTH 
C0NSTANT. 

ARRAY  NUMBER,  XXXXX,  IS  THE  DUPLICATE  0F  THE  XXXXXTH  ARRAY. 
FIXED  C0NSTANT  NAME  N0T  IN  LIST. 

NUMBER  0F  GEN  ARGUMENTS,  XXX,  EXCEEDS  NUMBER  REQUIRED. 

ST0RAGE  ALL0TTED  F0R  THIS  DATA  BL0CK  HAS  BEEN  EXCEEDED. 
PR0CESSING  WILL  RESUME  WITH  THE  NEXT  DATA  BL0CK. 

Forming  PCS— 

N0DE,  XXXXX,  HAS  N0  MATCH  IN  THE  NA-NB  PAIRS. 

ADJ0INING  NODE,  XXXXX,  0F  NA-NB  PAIR  HAS  N0  MATCH  IN  THE  N0DAL 
BL0CK,  C0NDUCTOR  IS  N0.,  XXXXX 

Processing  Program  Blochs— 

EXECUTI0N  BL0CKS  IN  IMPR0PER  0RDER  0R  ILLEGAL  BL0CK 
DESIGNATI0N  ENC0UNTERED. 

Explanation:  Some  alpha  character  other  than  K  or  A  has  been  used  to  reference  a 
data  block.  In  a  thermal  problem  a  designator  other  than  G,  K,  or  A 
is  assumed  to  be  referencing  the  nodal  block. 

MISSING  N0DE  NUMBER,  XXXXX. 

MISSING  C0NDUCT0R  NUMBER,  XXXXX. 

MISSING  C0NSTANT  NUMBER,  XXXXX. 

MISSING  ARRAY  NUMBER,  XXXXX. 

FIXED  C0NSTANT  NAME,  AAAAA,  N0T  IN  LIST. 
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NUMBER  0F  SUBR0UT1NES  REQUESTED  EXCEEDS  75. 
Explanation:  More  than  75  unique  subroutines  have  been  called. 


Processing  Parameter  Changes— The  first  five  parameter  change  error  messages  are 
prefaced  with  the  words:  PARAMETER  CHANGE  ERR0R. 

N0DE  NUMBER,  XXXXX,  WAS  N0T  DEFINED  IN  THE  0RIGINAL  PR0BLEM. 

C0NDUCT0R  NUMBER,  XXXXX,  WAS  N0T  DEFINED  IN  THE  0RIGINAL 
PR0BLEM. 

C0NSTANT  NUMBER,  XXXXX,  WAS  N0T  DEFINED  IN  THE  0RIGINAL 
PR0BLEM. 

ARRAY  NUMBER,  XXXXX,  WAS  N0T  DEFINED  IN  THE  0RIGINAL  PR0BLEM. 

C0NSTANTS  BL0CK  WAS  EMPTY  IN  THE  0RIGINAL  PR0BLEM. 

ARRAY  BL0CK  WAS  EMPTY  IN  THE  0RIGINAL  PR0BLEM. 

ARRAY  NUMBER  XXXXX  -  DIMENSI0NS  N0T  EQUAL.  0RIGINAL, 

XXXXX,  CHANGE,  XXXXX. 


Terminations  Due  to  Errors  (no  preceding  asterisks)— 

THE  AB0VE  PARAMETER  CHANGE  WILL  N0T  BE  EXECUTED. 

ERR0R  TERMINATI0N  -  L0ADING  IS  SUPPRESSED. 

V.  OPERATING  SYSTEM  DESCRIPTION 
General 

The  CDC-3800  (FORTRAN  IV)  version  of  CINDA  exists  logically  as  a  preprocessor, 
processor,  and  library.  The  operational  continuity  of  these  portions  is  made  possible  by 
the  CDC  Drum  SCOPE  system  (see  Fig.  5). 

The  function  of  the  preprocessor  is  to  operate  on  a  user-supplied  problem  and  pro¬ 
duce  the  following  items. 

1.  Processor  Main  Program— This  small  routine  acts  primarily  as  a  communications 
link  in  providing  addressing  relationships  between  the  operational  user  program  and  user 
data. 


2.  User  Program— These  FORTRAN  source  subroutines  are  operational  equivalents 
of  the  user’  Execution,  Variables  1  and  2,  and  Output  Calls  blocks. 

3.  User  Data— Binary  data  generated  consists  of  definitions  of  parameters  referenced 
in  the  various  user  data  blocks  and  their  corresponding  values. 
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Fig.  5— Flow  of  C1NDA  operating  system 

34 


NRt,  REPORT  7656 


Tiie  preprocessor  and  appropriate  use  of  the  CDC-3800  system  control  cards  allows  con¬ 
struction  of  the  above  from  tape  when  the  RECALL  option  is  utilized. 

The  processor  performs  reading  of  the  user  data  values  prepared  previously  and  calls 
the  user  program  (i.e.,  Execution  block). 

The  C1NOA  library  contains  a  large  number  of  various  types  of  subprograms  to  ac¬ 
complish  most  user  requirements.  Drum  SCOPE’S  LIBED1T  provides  simple,  flexible 
methods  for  the  maintenance  of  this  library.  In  addition,  it  is  not  necessary  that  a  sub¬ 
routine  be  updated  to  the  library  prior  to  availability  in  the  user  problem. 


Preprocessor 

Operation  of  the  Preprocessing  Phase.— (See  Fig.  6  for  flowchart.)  The  main  pro¬ 
gram  PREPRO  accomplishes  the  initialization  of  data  values  and  tape  units  and  defines 
the  order  of  processing  by  calling  seven  subroutines. 

1.  If  the  problem  being  processed  is  a  RECALL  problem,  subroutine  SPLIT  is 
called  to  read  the  recalled  problem  data  and  number  definitions  from  the  input  tape  and 
write  these  on  the  appropriate  work  tape.  SPLIT  calls  SKIP  if  the  input  tape  is  not  posi¬ 
tioned  at  the  problem  being  recalled  (see  Store  and  Recall  Problems  Options). 

2.  CODERD  reads  the  title  block  and  the  block  title  cards.  It  then  calls  DATARD, 
which  reads  the  free-form  data  cards  in  the  four  (or  two,  if  General  Problem)  data  blocks 
and  2ny  parameter  change  data.  Each  card  is  read,  a  format  is  constructed  for  it,  and 
then  it  is  reread.  The  data  from  each  block  are  written  on  the  data  tape  as  one  record. 
The  number  definitions  of  the  data  and  the  NA-NB  pairs  are  written  on  work  tapes. 

3.  PSEUDO  reads  the  node  number  definitions  and  NA-NB  pairs  from  work  tape. 
The  PCS  (long  or  short)  is  constructed,  packed  by  PACK43  and  flagged  by  ORMIN,  and 
written  on  the  data  tape.  PACK43  and  ORMIN  call  BIT,  a  COMPASS  packing  and  un¬ 
packing  routine. 

4.  GENLNK  constructs  the  main  program  of  the  processor  (LINKO),  including 
COMMON  and  DIMENSION  information.  BLKCRD  and  STFFB  are  called  upon  to  fill 
an  array  with  FORTRAN  source  code,  which  is  then  written  onto  logical  unit  LB4P  by 
WRTBLK.  WRTSCOPE  writes  SCOPE  at  the  end  of  LB4P  after  completion  of  the  other 
four  subroutines  (EXECTN,  VARBL1,  VARBL2,  and  OUTCAL). 

5.  PRESUB  reads  the  title  cards  of  the  four  program  blocks  and  initiates  the  con¬ 
struction  of  each  new  subroutine.  CINDA4  converts  the  CINDA  "calls”  in  the  program 
blocks  into  FORTRAN  subroutine  calls.  Data  referenced  by  input  number  definition  is 
changed  to  refer  to  its  relative  location  in  COMMON  data  arrays. 

6.  INITAL  combines  the  original  set  of  the  data  and  the  nitial  parameter  changes 
and  writes  the  updated  set  of  data  on  the  data  tape. 
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7.  FINAL  converts  final  parameter  change  data  (number  definitions  and  values)  to 
relative  array  locations  and  values  and  writes  number-value  records  on  the  data  tape. 
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VI.  EXTERNAL  PLOTTING  PACKAGE 

CINDA’s  plot  package,  for  use  on  the  CALCOMP  plotter,  is  an  external  program  that 
will  plot  a  graph  of  time  vs  temperature  for  each  problem  node.  The  input  to  *his  pro¬ 
gram  is  an  output  file  (unit  24)  generated  by  TSAVE  during  a  previous  CINDA  problem 
run.  The  program  can  be  run  separately  from  the  CINDA  problem  using  the  tape  from 
TSAVE  as  input,  or  it  can  be  placed  behind  the  CINDA  problem  in  a  single  run.  In  the 
latter  case,  unit  24  may  be  either  a  drum  unit  or  a  tape  equipped  for  later  use. 

The  package,  available  as  a  binary  deck,  consists  of  three  routines.  The  main  pro¬ 
gram,  PLOTTEMP,  calls  PLOTPREP,  which  rearranges  the  data  from  the  input  tape  and 
writes  it  on  unit  25.  Unit  25  contains  the  actual  node  numbers,  the  time  array,  and  the 
temperature  profile  for  each  node.  PLOTTEMP  then  reads  a  set  of  data  cards  which  give 
the  plot  heading,  X-  and  Y-axis  limits,  and  nodes  to  be  plotted.  The  temperature  array 
for  each  node  is  read,  and  if  the  node  is  to  be  plotted,  PLOTT  is  called,  which  in  turn 
activates  CALCOMP  routines.  A  separate  set  of  axes  is  drawn  for  each  node.  When  all 
temperatures  have  been  read,  the  tape  is  rewound  and  a  new  set  of  data  (if  any)  is  read. 


Data  Set 

A  data  set  consists  of  at  least  three  cards. 

♦CARD  1 -TITLE 

Columns  1—40  will  be  used  as  the  plot  heading. 

♦CARD  2-AXIS  LIMITS  and  TEMPERATURE  SCALE  OPTION 

Four  floating-point  values  musl  be  entered  in  an  E9.2  format.  The  minimum 
and  maximum  times  (X-axis)  must  be  in  fields  4-12  and  14-22,  respectively;  the  mini¬ 
mum  and  maximum  temperatures  (Y-axis)  must  be  in  fields  24-32  and  34-42,  respectively. 
Column  43  contains  the  temperature  scale  option.  A  blank  indicates  that  the  data  will  be 
p'otted  in  Fahrenheit  temperatures,  and  a  1  specifies  Centigrade. 

4  12  22  32  43 

fill  1 

NNNNNN.NN  NNNNNN.NN  NNNNNN.NN  NNNNNN.NN  I 

(TIM EMIN)  (TIMEMAX)  (TEMPMIN)  (TEMPMAX)  (Centigrade  scale) 

♦CARD  3  -NODES  and  OPTION 


Options: 

Plot  all  nodes— card  is  blank. 

Plot  certain  nodes— column  1  contains  a  nodes  are  listed. 
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Plot  all  bat  certain  nodes— column  1  contains  any  character  but  *  or  $; 
nodes  are  listed. 

The  nodes  may  be  listed  individually  or  in  inclusive  pairs.  Each  node  must  be  fol¬ 
lowed  by  a  comma  except  for  the  inclusive  pair,  which  must  be  separated  by  a  blank  in¬ 
stead  of  a  comma.  (The  second  node  of  the  pair  must  be  followed  by  a  comma,  how¬ 
ever.)  An  asterisk  (*),  instead  of  a  comma,  must  follow  the  last  node. 

The  node  numbers  must  be  in  14  format,  right  adjusted  to  columns  5,  10,  15,  .  .  ., 
80.  The  commas  or  separating  blanks  go  in  columns  6,  11,  16,  .  .  .,  76.  Data  may  be 
continued  to  a  following  card  by  putting  a  comma  or  blank  in  column  1  of  that  card. 
Data  cards  will  be  read  until  an  asterisk  is  encountered  (limit  of  nine  cards). 


Examples 


1.  1 

5 

10 

15 

20 

1 

1 

1 

1 

1 

$ 

10, 

15, 

25 

30* 

Plot  nodes  10,  15,  25,  26,  27,  28,  29,  30.  (Pairs  must  be  in  ascending  order; 
otherwise,  order  of  magnitude  isn’t  important.) 


to 

5 

10 

15 

20 

25 

1 

1 

1 

1 

1 

1 

+ 

5, 

8, 

10 

12, 

15* 

Plot  all  nodes  but  5,  8,  10,  11,  12,  15. 


3.  1  5  10  15  20  .  75  80 

11111  11 

(Card  1)  $1050,1060  1062,1070 .  2000,2005 

(Card  2)  2007,2012* 


Plot  1C50, 1060,1061,1062,1070, ...  ,2000,2005,2006,2007,2012. 
(A  comma  in  column  1  of  card  2  would  exclude  node  2006.) 


4.  1  5  10  15  80 

1111  1 

(Card  1)  +  30,  25  26, .  100 

(Card  2)  * 


Plot  all  but  nodes  30,  25,  26,  27,  28,  ....  100. 

Another  set  of  data  may  follow  the  last  node  card,  thus  enabling  the  programmer  to 
redefine  the  title  and  limits  and  to  plot  different  nodes.  Any  number  of  data  sets  may 
be  used. 
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Diagnostics 

Messages  that  may  be  encountered  while  plotting  are  listed  below.  The  first  two  will 
not  cause  job  termination. 

1.  ***  S0ME  NODES  T0  BE  PLOTTED  WERE  N0T  F0UND  0N  INPUT  TAPE, 
0R  WERE  DUPLICATED  0N  N0DE  CARD.  *** 

N0DES  0N  INPUT  TAPE  - 
(listing  of  nodes) 

N0DES  T0  BE  PLOTTED— 

(listing  of  nodes) 

Cause:  The  number  of  nodes  to  be  plotted  is  larger  than  the  number  of  nodes 
actually  plotted. 

2.  IN  F0LL0WING  N0DE  CARD,  AN  INCLUSIVE  PAIR  0F  N0DES  IS 
FOLLOWED  BY  A  BLANK-A  C0MMA  IS  ASSUMED  AND  PROCESSING 
IS  CONTINUED. 

3  ******  wrong  FORMAT  USED  0N  FOLLOWING  NODE  CARD- 

PROCESS  NEXT  SET  OF  DATA,  IF  ANY. 

Cause:  A  character  other  than  a  comma,  blank,  or  asterisk  has  been  found  in  a 
column  intended  for  those  characters  only. 

This  error  terminates  the  processing  of  a  current  data  set,  and  processing  con¬ 
tinues  to  the  next  set  if  there  is  one.  (Most  data  format  errors  will  cause  job 
termination  by  the  system.) 


VII.  TAPE  USAGE  AND  DECK  SETUPS 

This  section  shows  deck  setups  and  t3pe  usage  for  various  types  of  runs  on  the  3800 
Drum  SCOPE  system. 


CINDA  Operating  System 

Table  1  lists  the  program  and  data  files  used  in  the  preprocessor,  processor,  and 
library. 

Units  15,  17,  18,  19,  and  27  are  preprocessing  units  only  and  are  available  as  scratch 
units  during  processing.  Units  16,  21,  22,  24,  30,  33,  and  40  can  also  be  used  for  that 
purpose  if  the  corresponding  options  are  not  activated. 
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Table  1 

CINDA  Tape  Usage 


Logical 

Unit 

Program 

Variable 

Function 

*9 

— 

CINDA  Master  tape;  (file  1  contains  preproc¬ 
essor,  file  2  contains  Library). 

12 

LB3D 

Data  tape  (original  problem  and  all  parameter 
changes). 

14 

LB4P 

Program  tape  (contains  generated  FORTRAN 
routines  LINK0,  EXECTN,  VARBL1, 

VARBL2,  OUTCAL). 

15 

LUT7 

Variables  1  calls  generated  from  node  and 
conductor  data  blocks. 

**16 

— 

Matrix  retrieval  unit  in  library. 

17 

NA-NB  pairs;  data  number  definitions. 

(From  parameter  changes.) 

18 

LUT3 

Copy  of  original  problem  data. 

19 

LUT4 

Parameter  Change  data. 

20 

LUT1 

Data  number  definitions. 

**21 

— 

Problem  recall  data  tape. 

**22 

— 

Problem  store  data  tape. 

**24 

— 

Output  from  TSAVE. 

27 

INTERN 

Data  conversion  scratch  tape. 

**30 

KRR 

FORTRAN  reread  unit  in  preprocessor. 

Matrix  storage  unit  in  library. 

33 

— 

Scratch  unit  in  STOREP. 

**40 

Binary  program  tape  (processor)  used  with 
store  and  recall  options. 

♦Equipped  units. 

•♦Equipped  units  depending  on  options.  Matrix  storage  and  retrieval  requires  equipping  tapes  16  and  30. 
The  STOREP  option  requires  equipping  tapes  40  and  22;  the  RECALL  option  requires  40  and  21  (and 
22  if  desired). 


General  deck  structures  for  different  kinds  of  CINDA  runs  are  shown  below.  The 
character  A  denotes  a  7  and  9  punch  in  the  column,  and  a  A  is  for  an  ll(-),  0,  7,  9 
punch.  The  number  enclosed  in  parentheses  in  the  job  card  (e.g.,  5)  may  have  to  be  in¬ 
creased  if  several  tapes  are  used.  The  current  label  for  tape  unit  9  is  CINDA  MASTER, 
1,1,999.  The  other  tapes  may  be  unlabeled.  Corresponding  job  request  forms  are  found 
in  Fig.  7. 
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Not  RECALL  (Ordinary  Run)— (See  Fig.  8.) 
AJ0B(5), ... 

AEQUIP,  9=  (label),  R0,  HI,  DA 
ArBANK,(0),/4/ 

AL0AD.9 
ARUN.t.S 
blank  card 


problem  data  deck  through  BCD  3END  0F  DATA 


ALIBR  AR  Y,  9,  LCIND  A 
AFTN,  I=14,L,X 
AL0AD 

binary  decks,  if  any 


• 

ARUN,t,G 

E0F 

Not  RECALL  (STOREP  run) 
AJ0B(5), ... 

AEQUIP,  9=  (label),  R0,  HI,  DA 
AEQUIP,40= (label),  RV.',  HI,  DA 
AEQUIP, 22=  ( label), W0, HI, DA 
ABANK,(0),/4/ 

AL0AD.9 
ARUN,t,  1 
blank  card 


problem  data  deck  through  BCD  3END  0F  DATA  (includes  at  least  one  call  to 
STOREP) 


ALIBRARY,9,LCINDA 

AFTN,I=14,L,X=40 

AL0AD.4O 


binary  decks,  if  any 


ARUN,t,e 

E0F 
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RECALL 

AJ0B(5), . . . 

AEQUIP,9=(label),R0,HI,DA 

AEQUIP,4O=(label),R0HI,DA 

AEQUIP,21=(label),R0,HI,DA 

*AEQUIP,22=(label),YV0,HI,DA 

ABANK,(0),/4/ 

AL0AD.9 
ARUN,t,C 
RECALL  Card 


INITIAL  PARAMETERS  blocks  and  BCD  3END  0F  DATA 


ALIBRARY,9,LCINDA 

AL0AD.4O 


binary  decks  if  any 


ARUN,t,fi 

E0F 

For  any  other  options  using  tapes,  the  tapes  should  be  equipped  as  those  shown 
above,  using  the  appropriate  logical  unit  number  and  label.  The  user  must  also  designate 
whether  the  tape  is  read  only  (R0),  only  written  on,  (W0),  or  both  (RW).  In  the  latter 
case,  if  the  tape  is  written  on  first,  '  .e  output  block  of  the  job  request  form  is  checked. 
If  the  tape  is  read,  then  written  on,  both  the  input  and  output  blocks  should  be  checked. 
(Check  the  Drum  SCOPE  manual  for  more  details.) 


Plot  Package 

Table  2  lists  the  files  used  in  the  plot  packages. 


‘Optional— used  only  if  problem  is  to  be  restored. 
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Table  2 


Logical  Unit 

Function 

10 

Plotting  unit 

24 

Output  from  TSAVE; 

Input  to  PLOTPREP 

25 

Output  from  PLOTPREP; 

Input  to  PLOTT 

NOTE:  Units  25  and  10  are  drum  units.  Unit  24  is  usually  a  tape, 
but  can  also  be  a  drum  unit  if  the  plotting  run  follows  the  CINDA 
run  in  the  same  job. 


Below  are  sample  deck  sets  showing  a  CINDA  run  that  generates  the  TSAVE  tape, 
and  a  plotting  run  that  uses  the  TSAVE  tape  as  input.  Figure  9  shows  a  job  request 
form  for  the  plotting  run.  Figure  10  displays  the  deck  setup  for  a  combined  CINDA  and 
plotting  run. 

Not  RECALL  (Generate  TSAVE  tape) 

AJ0B(5), . . . 

AEQUIP,9=(label),R0,HI,DA 

AEQUIP,24=(label),YV0,HI,DA 

ABANK,(0),/4/ 

AL0AD.9 
ARUN,t,J! 
blank  card 


problem  data  deck  through  BCD  3END  0F  DATA  (includes  a  call  to  TSAVE 
in  Output  Calls) 


ALIBRARY,9,LCINDA 

AFTN,1=14,L,X 

AL0AD 


binary  decks,  if  any 


ARUN,t,e 

EOF 


MARY  E.  GEALY 


Fig.  9— Sample  Job  Request  form  for  plot  run 


Plotting  Run 


AJ0B, .... 

AEQUIP,24=(label),R0,HI,DA 


binary  deck-< 


PL0TTEMP 

PL0TPREP 

PL0TT 


ARUN,t,l 


plotting  data 
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VIII.  ALPHABETIC  LISTING  OF  AVAILABLE  SUBROUTINES 


Name 

Page 

Name 

Page 

Name 

Page 

5 

j 

AABB 

CMPYI 

D11MDA 

k 

£ 

ABLATS 

CNBACK 

D11MDI 

* 

% 

ACSARY 

CNEXPN 

D12CYL 

-i 

*5 

ADARIN 

CNFAST 

D12MCY 

ADD 

CNFRWD 

D12MDA 

X 

X 

ADDALP 

CNFVVBK 

D2DEG1 

i 

f 

ADDARY 

C0LMAX 

D2DEG2 

l 

ADDFIX 

C0LMIN 

D2D1WM 

l 

ADDINV 

C0LMLT 

D2D2WM 

ALPHAA 

C0PY 

D2MXD1 

r, 

\ 

ARCC0S 

C0SARY 

D2MXD2 

ARCSIN 

CSGDMP 

D2MX1M 

ARCTAN 

CSQR.I 

D2MX2M 

\ 

ARINDV 

CVQ1HT 

D3DEG1 

<. 

ARYADD 

CVQ1WM 

D3D.1WM 

ARYDIV 

DA11CY 

EFACS 

ARYEXP 

DA11MC 

EFASN 

l 

ARYINV 

DA12CY 

EFATN 

i 

ARYMNS 

DA12MC 

EFC0S 

: 

ARYMPY 

DIAG 

EFEXP 

\ 

ARYPLS 

DISAS 

EFFEMS 

ARYST0 

DIVARY 

EFFG 

ARYSUB 

DIVFIX 

EFL0G 

ASNARY 

DIVIDE 

EFP0VV 

ASSMBL 

D1DEG1 

EFSIN 

ATNARY 

D1DEG2 

EFSQR. 

' 

BIT 

D1DG1I 

EFTAN 

i 

BIVLV 

D1D1DA 

ELEADD 

* 

BKARAD 

D1D1IM 

ELEDIV 

BLDARY 

D1D1MI 

ELEINV 

* 

BRKARY 

D1D1WM 

ELEMUL 

BTAB 

D1D2DA 

ELESUB 

i 

BVSPDA 

D1D2WM 

ENDM0P 

BVSPSA 

D1MDG1 

E0F 

■i 

CALL 

D1MDG2 

EXPARY 

CDIVI 

D1M1DA 

EXPNTL 

t 

CINC0S 

D1M1MD 

FILE 

CINDSL 

D.1M1WM 

FIX 

CINDSM 

D1M2DA 

FLIP 

CINDSS 

D1M2MD 

FL0AT 

i 

CINSIN 

D1M2WM 

GENALP 

CINTAN 

D11CYI 

GENARY 

CMPXDV 

D11DAI 

GENC0L 

: 

CMPXMP 

D11DIM 

GENM 

CMPXSR 

D11MCY 

GENST 

48 


n 

- > ^  4  . .  11  l 

% 

k 

ft 

\ 

£ 

| 

£ 

s- 

Tt 

? 

3 

p 
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4 

p 

a 

? 

£ 

Name  Page 

Name  Page 

Name  Page  .■ 

i* 

»: 

c 

GSL0PE 

PRINT 

SMPINT  ] 

s 

A 

HEDC0L 

PRINTA 

SPLIT  | 

f 

INPUTG 

PRINTL 

SPREAD  | 

SPRESS  1 

INPUTT 

PRNTMA 

1 

INTRFC 

PRNTMP 

SQR00T  % 

1NVRSE 

PSINTR 

SQR0TI  | 

y 

# 

IRRADE 

PSNTWM 

STATE 

?. 

iRRADI 

PS0FTS 

STFSEP  1 

« 

ITRATE 

PUNCH 

STFSEQ  \ 

r 

JAC0BI 

PUNCHA 

STFSQS  S 

j* 

JC  IN 

PYMLT1 

STIFF  \ 

it 

LAGRAN 

QF0RCE 

STNDRD  t 

S 

j 

LGRNDA 

QINTEG 

ST0ARY  \ 

5c 

* 

LINE 

QINTGI 

ST0REP  I 

7 

$ 

LIST 

QMAP 

SUB  ! 

jr 

r< 

L0GE 

QMETER 

SUBARY  1 

1 

rJ 

L0GEAR 

QMTRI 

SUBFIX  \ 

3 

> 

L0GT 

RDTNQS 

SUMARY  j 

SYMINV  5 

£ 

f 

L0GTAR 

READ 

t 

LQDVAP 

REFLCT 

SYMLST  1 

t 

LSTAPE 

REWIND 

TANARY  ) 

3 

LSTPCS 

R0WMLT 

T0PLIN  i 

| 

LSTSQU 

RTP0LY 

TPRINT  \ 

ft 

MASS 

SCALAR 

TRANS  t 

e 

MATRIX 

SCALE 

TRNPRT  ? 

MAXDAR 

SCLDEP 

TRPZDA  l 

| 

MLTPLY 

SCLIND 

TRPZD  1 

M0DES 

SCRPFA 

TSAVE  ? 

l 

MPYARY 

SETMNS 

TS0FP  ? 

u 

MPYFIX 

SETPLS 

UNPAK  1 

•- 

i* 

MULT 

SETUP 

UPDM0P  ^ 

bj 

MXDRAL 

SHFTV 

UNITY  i 

V 

NEWRT4 

SHFTVR 

VARCCM  1 

} 

NEWTRT 

SHIFT 

VARCSM  1 

<5. 

0NES 

SHUFL 

VARC1  I 

1 

PLYARY 

SIGM* 

VARC2  1 

i 

* 

PLYEVL 

SIMEQN 

VARGCM  | 

PLYNML 

SINARY 

VARGSM  I 

PNTABL 

SKPLIN 

VARG1  | 

j’’ 

P0LMLT 

SLDARD 

VARG2  1 

U 

P0LS0V 

SLDARY 

WRITE  1 

? 

fc 

P0LVAL 

SLRADE 

WRTARY  I 

5 

P0LYADD 

SLRADI 

WRTL08  1 

f 

? 

< 

5 

t 

6 

PRESS 

SM0PAS 

ZER0  | 

t, 

- 

> 

49 

Sf 

V 

■ss 

Jti 

* 

? 

f, 

*5 

tr 

X 

j; 

5 

i 

! 

If 

J 

7 

■L. 

.  Jj 

Execution  Subroutines 
Name 


CINDSS  (Steady  stole 

C1NDSL  (Steady  state 

CINDSM  (Steady  state 

CNFRWD  (Explicit  forward  differencing) 

CNFAST  (Accelerated  forward  diffe^cififl^j^,, t 

CNEXPN  (Explicit  exponential  prjwi.£),ion  f  %f 

CNFWBK  (Implicit  forward -back  wtc^^riferencitig) 

CNBACK  (Implicit  backward  differencing) 


Execution  Subroutine  CINDSS 


-<'i 

• 

Purpose— This  subroutine  ignores  the  capacitance  valuer  &&g?P<ra&iV?  Hu#*-  tv’  . 
late  the  network  steady  state  solution.  Due  to  the  SPOS  re^ibw-nsrpi', 
solved  by  a  “block”  iterative  method.  However,  if  ail  diffuser  tukitfS  japsre  toSS 

arithmetic  nodes  they  would  be  calculated  by  a  successive  j «£ffeUv 


**  | 

.iipK  < 
v  l  A. 

AV  '  ;‘l1 

'•  '-h^aa- 

>  •*  ft  u 


ARLXCA  for  arithmetic  nodes).  The  subroutine  will  comv.a ;i)V:-  v-,j  iwhn$»  \  r  >T  r  ’•fef6  't. 

above  criteria  is  met.  If  the  iteration  count  exceeds  NLOO-'-to;?  ftjspw$n«.&1 
printed.  Variables  1  and  Output  Calls  are  performed  ..the  sfetp} 

put  Calls  are  performed  upon  completion.  If  not  specified,  cvuVol  •  V  -  '$3 

DAMPA  are  set  at  1.0.  They  are  used  as  multipliers  times  tf«,  ‘  ‘  <v'*  '  1 

1.0  minus  their  value  is  used  as  multipliers  times  the  old  tenq 
the  returned  answer.  This  weighting  of  so  much  new  ftp*!  so 
ing  oscillations  due  to  nonlinearities.  Thc>  mav  also  he 

If  a  series  of  steady  state  solutions  at  various  times 
by  specifying  control  constants  TIM  END 
the  output  interval  and  the  computation 
have  to  be  made  in  Variables  1  to  modify 

If  desired,  the  CINDSS  call  can  he  followed 
subroutines  which  has  the  same  SPCS  requiremirsiS) 

tion  becomes  the  initial  conditions  for  the  transient)  fKUwfvet  fiKl* 

utilizes  control  constants  TIMEND  and  OUTPUT  the  tjsez, m,g. ,sjy?edy  Wic* S^&- 
execution  block  after  the  steady  state  call  and  pnor  to  ike  vjmu&ii  spiifyti.x  -yi-  *,.• 


Restrictions— The  SPCS  option  is  required..  Diffusion  node?  receive  a 
tion,  while  arithmetic  nodes  receive  *»***'*»•* 

are  utilized.  Control  constants  NLOOP 
fied.  Successive  steady  state  solutions 
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TIMEND  and  OUTPUT.  Other  control  constant;  which  are  activated  or  used  are 
LOOPCT,  DRLXCC  and/or  ARLXCC,  TIMEN,  TIMEM,  T1MEO,  DAM  PD,  DAMPA, 
DTIMEU,  L1NECT,  and  PAGECT.  Control  constant  OPE1TR  is  checked  for  output  each 
iteration. 

Calling  Sequence— C INDSS-  This  subroutine  utilizes  one  dynamic  storage  core  loca¬ 
tion  for  each  diffusion  node. 


Execution  Subroutine  CINDSL 

Purpose— This  subroutine  ignores  the  capacitance  values  of  diffusion  nodes  to  calcu¬ 
late  the  network  steady  state  solution.  Since  this  subroutine  has  the  LPCS  requirement, 
both  diffusion  and  arithmetic  nodes  receive  a  successive  point  iteration.  In  addition,  each 
third  iteration  a  linear  extrapolation  is  performed  on  the  error  function  plot  of  each  node 
in  an  attempt  to  accelerate  convergence.  The  user  is  required  to  specify  the  maximum 
if  number  of  iterations  to  be  performed  in  attempting  to  reach  the  steady  state  solution 
-  (control  constant  NLOOP)  and  the  relaxation  criterion,  which  determines  when  it  has 
been  reached  (DRLXCA  for  diffusion  nodes  and/or  ARLXCA  for  arithmetic  nodes).  The 
subroutine  will  continue  to  iterate  until  one  of  the  above  criteria  is  met.  If  the  iteration 
count  exceeds  NLOOP  an  appropriate  message  is  printed.  Variables  1  and  Output  Calls 
are  performed  at  the  start,  and  Variables  2  and  Output  Calls  are  performed  upon  comple¬ 
tion.  If  not  specified,  control  constants  DAMPD  and  DAMPA  are  set  at  1.0.  They  are 
used  as  multipliers  times  the  new  temperatures  while  1.0  minus  their  value  is  used  as 
multipliers  times  the  old  temperatures  in  order  to  weight  the  returned  answer.  This 
weighting  of  so  much  new  and  so  much  old  is  useful  for  damping  oscillations  due  to  non- 
lir.earities.  They  may  also  be  used  to  achieve  overrelaxation. 

If  a  series  of  steady  state  solutions  at  various  times  is  desired  it  can  be  accomplished 
by  specifying  control  constants  TIMEND  and  OUTPUT.  OUTPUT  will  be  used  both  as 
the  output  interval  and  the  computation  interval.  In  this  case  appropriate  calls  would 
have  to  be  made  in  Variables  1  to  modify  boundary  conditions  with  time. 

If  desired,  the  CINDSL  call  can  be  followed  by  a  call  to  one  of  the  transient  solution 
subroutines  which  has  the  same  LPCS  requirement.  In  this  manner  the  steady  state  solu¬ 
tion  becomes  the  initial  conditions  for  the  transient  analysis.  However,  since  CINDSL 
utilizes  control  constants  TIMEND  and  OUTPUT  the  user  must  specify  their  values  in  the 
execution  block  after  the  steady  state  call  and  prior  to  the  transient  analysis  call. 


Restrictions— The  LPCS  option  is  required.  Diffusion  and  arithmetic  nodes  receive  a 
successive  point  iteration  and  an  extrapolation  method  of  acceleration.  Control  constants 
NLOOP  and  DRLXCA  and/or  ARLXCA  must  be  specified.  Successive  steady  state 
solutions  can  be  obtained  by  specifying  control  constants  TIMEND  and  OUTPUT.  Other 
control  constants  which  are  activated  or  used  are:  LOOPCT,  DRLXCC,  and/or  ARLXCC, 
TIMEN,  TIMEM,  T1MEO,  DAMPD,  DAMPA,  DTIMEU,  LINECT,  and  PAGECT.  Control 
constant  OPE1TR  is  checked  for  output  each  iteration. 
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Calling  Sequence  -CINDSL  —  This  subroutine  utilizes  two  dynamic  storage  core  loca¬ 
tions  for  each  diffusion  and  arithmetic  node. 


Execution  Subroutine  CINDSM 

Purpose— This  subroutine  is  designed  to  calculate  the  network  steady  state  solution 
of  moderately  radiation-dominated  problems.  It  is  similar  to  CINDSL  in  that  the  LPCS 
option  is  required  and  that  all  nodes  receive  a  successive  .point  iteration  and  the  same  ex¬ 
trapolation  method  of  acceleration.  Other  execution  subroutines  evaluate  the  nonlinear 
radiation  conductors  each  time  they  are  encountered  during  an  iteration.  CINDSR  differs 
in  that  it  linearizes  the  nrcblem  by  calculating  effective  radiation  conductors  and  solves 
the  linearized  problem.  It  then  reevaluates  the  effective  radiation  conductors,  solves  the 
linear  problem  and  continuously  repeats  the  process.  The  user  must  specify  the  maximum 
number  of  iterations  to  perform  in  attempting  to  reach  the  steady  state  solution  and  the 
energy  balance  of  the  system  to  be  satisfied  as  a  criterion.  This  system  energy  balance  is 
the  difference  between  all  energy  into  the  system  and  all  energy  out  and  is  specified  as 
control  constant  BALENG.  CINDSM  internally  calculates  the  iterative  relaxation  criteria 
damping  factors  and  loopings  to  be  performed  in  solving  the  linearized  problem.  It  con¬ 
tinuously  increases  the  severity  of  the  relaxation  criteria  until  the  BALENG  criteria  is 
met  for  two  successive  linearized  problems  with  virtually  no  temperature  change  between 
the  two.  Systems  with  small  energy  transfer  rates  to  the  boundaries  are  difficult  to  solve. 
A  reasonable  rule  is  to  set  BALENG  at  1 %  of  the  rate  in  or  out.  Successive  steady  state 
analyses  may  be  performed  and  CINDSM  may  be  followed  by  a  call  to  a  transient  analysis 
routine  with  the  same  LPCS  option  requirement. 

Restrictions— The  LPCS  option  is  required.  Control  constants  NLOOP,  LAXFAC, 
and  BALENG  must  be  specified  and  be  greater  than  zero.  DAMPD  may  be  used.  If  it  is 
not  specified,  the  routine  will  set  DAMPD  to  1.0.  Successive  steady  state  solutions  can 
be  obtained  by  specifying  control  constants  TIMEND  and  OUTPUT.  Other  control  con¬ 
stants  which  are  activated  or  used  are  LOOPCT,  ENGBAL,  and/or  ARLXCC,  TIMEN, 
TIMEM,  TIMEO,  DTIMEU,  LINECT,  and  PAGECT.  Control  constant  OPE1TR  is  checked 
for  output  each  iteration.  Caution:  Each  radiation  conductor  must  have  a  unique  con¬ 
ductor  number. 


Calling  Sequence— CINDSM  —  This  subroutine  utilizes  three  dynamic  storage  core 
locations  fo *•  each  diffusion  and  arithmetic  node  and  one  more  for  each  radiation  cor 
conductor. 


Execution  Subroutine  CNFRWD 

Purpose— This  subroutine  performs  transient  thermal  analysis  by  the  explicit  forward- 
differencing  method.  The  stability  criterion  of  each  diffusion  node  is  calculated  and  the 
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minimum  value  is  placed  in  control  constant  CSGMIN.  The  time  step  used  (control  con¬ 
stant  DTIMEU)  is  calculated  as  95%  of  CSGMIN  divided  by  CSGFAC.  Control  constant 
CSGFAC  is  set  at  1.0  unless  specified  larger  by  the  user.  A  “look-ahead”  feature  is  used 
when  calculating  DTIMEU.  If  one  time  step  will  pass  the  output  time  point,  the  time 
step  is  set  to  come  out  exactly  on  the  output  time  point;  if  two  time  steps  will  pass  the 
output  time  point,  the  time  step  is  set  so  that  two  time  steps  will  come  out  exactly  on 
the  output  time  point.  DTIMEU  is  also  compared  to  DTIMEH  and  DTIMEL.  If  DTIMEU 
exceeds  DTIMEH  it  is  set  equal  to  it,  if  DTIMEU  is  less  than  DTIMEL  the  problem  is 
terminated.  If  no  input  values  are  specified,  DTIMEL  is  set  at  zero  and  DTIMEH  it  is 
set  at  infinity.  The  maximum  temperature  change  calculated  over  an  iteration  is  placed 
in  control  constant  DTMPCC  and/or  ATMPCC.  They  are  compared  to  DTMPCA  and/or 
ATMPCA,  respectively,  and  if  larger  cause  DTIMEU  to  be  modified  so  that  they  com¬ 
pare  as  equal  to  or  less  than  DTMPCA  and/or  ATMPCA.  If  DTMPCA  and/or  ATMPCA 
are  not  specified  they  are  set  at  infinity. 

All  diffusion  nodes  are  calculated  prior  to  solving  the  arithmetic  nodes.  The  user 
may  iterate  the  arithmetic  node  solution  by  specifying  control  constants  NLOOP  and 
ARLXCA.  If  the  arithmetic  node  iteration  count  exceeds  NLOOP,  answers  are  ac¬ 
cepted  as  is  and  the  subroutine  continues  without  any  user  notificai.  >n.  In  addition  the 
user  may  specify  control  constant  DAMPA  in  order  to  dampen  possible  oscillations  due 
to  nonlinearities.  The  arithmetic  nodes  may  be  used  to  specify  an  incompressible  pressure 
or  radiosity  network.  In  this  manner  they  would  be  solved  implicitly  each  time  step,  but 
evaluation  of  temperature  varying  properties  would  suffer  a  lag  of  one  time  step. 


Restrictions— The  SPCS  option  is  required  and  control  constants  TIMEND  and 
OUTPUT  must  be  specified.  Problem  start  time  if  other  than  zero  may  be  specified  as 
TIMEO.  Other  control  constants  used  or  activated  are:  T1MEN,  TIMEM,  CSGMIN, 
CSGFAC,  DTIMEU,  DTIMEL,  DTIMEH,  DTMPCA,  DTMPCC,  ATMPCA,  ATMPCC, 
NLOOP,  LOOPCT,  DAMPA,  ARLXCA,  ARLXCC,  OPEITR,  BACKUP,  LINECT,  and 
PAGECT. 


Calling  Sequence— CNFRWD  —  This  subroutine  utilizes  one  dynamic  storage  core 
location  for  each  diffusion  and  arithmetic  node. 


Execution  Subroutine  CNFAST 

Purpose— This  subroutine  is  a  modified  version  of  CNFRWD  which  allows  the  user 
to  specify  the  minimum  time  step  to  be  taken.  The  time  step  calculations  proceed  ex¬ 
actly  as  in  CNFRWD  until  the  check  with  DTIMEL  is  made.  If  DTIMEU  is  less  than 
DTIMEL  it  is  set  equal  to  it.  As  each  node  is  calculated  its  CSGMIN  is  obtained  and 
compared  to  DTIMEU.  If  equal  to  or  greater,  the  nodal  calculation  is  identical  to 
CNFRWD.  If  the  CSGMIN  for  a  node  is  less  than  DTIMEU  the  node  receives  a  steady 
state  calculation.  If  only  a  small  portion  of  the  nodes  in  a  system  receive  the  steady 
state  calculation  the  answers  are  generally  reasonable.  However,  as  the  number  of  nodes 
receiving  steady  state  calculations  increases,  so  do  the  solution  inaccuracies. 
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Restrictions— The  SPCS  option  is  required  and  control  constants  TIMEND  and 
OUTPUT  must  be  specified.  The  checks  on  control  constants  DTMPCA,  ATMPCA  and 
BACKUP  are  not  performed.  Other  control  constants  which  are  used  or  activated  are 
T1MEN,  T1MEM,  TIMEO,  CSGMIN,  CSGFAC,  DTIMEU,  DTIMEL,  DTIMEH,  DTMPCC, 
ATMPCC,  DAMPA,  ARLXCA,  ARLXCC,  NLOOP,  LOOPCT,  LINECT,  and  PAGECT. 


Calling  Sequence— CNF  AST  —  This  subroutine  utilizes  one  dynamic  storage  core 
locatior  for  each  diffusion  node. 

Execution  Subroutine  CNEXPN 

Purpose— This  subroutine  performs  transient  thermal  analysis  by  the  exponential 
prediction  method,  and  the  solution  equation  is  of  the  following  form: 


This  equation  is  unconditionally  stable,  no  matter  what  size  time  step  is  taken,  and  it  re¬ 
duces  to  the  steady  state  equation  for  an  infinite  time  step.  However,  stability  is  not  to 
be  confused  with  accuracy.  Time  steps  larger  than  would  be  taken  with  CNFRWD  re¬ 
main  stable  but  tend  to  lose  or  gain  energy  in  the  system.  For  this  reason  this  subroutine 
is  not  recommended  where  accuracy  is  sought.  However,  it  is  suitable  for  parametric 
analysis  where  trends  are  sought  and  a  more  accurate  method  will  be  utilized  for  a  final 
analysis. 

The  inner  workings  of  the  subroutine  are  virtually  identical  to  CNFRWD  with  the 
exception  of  the  solution  equation  and  the  use  of  CSGFAC.  The  time  step  used 
(DTIMEU)  is  calculated  as  CSGMIN  times  CSGFAC.  The  look-ahead  feature  for  calcu¬ 
lating  the  time  step  is  identical,  as  are  the  checks  with  DTIMEH,  DTIMEL,  and  DTMPCA. 
The  diffusion  nodes  are  calculated  prior  to  the  arithmetic  nodes,  and  the  arithmetic  nodes 
utilize  NLOOP,  ARLXCA,  and  DAMPA,  exactly  the  same  as  CNFRWD. 

Restrictions— The  SPCS  option  is  required  and  control  constants  TIMEND  and 
OUTPUT  must  be  specified.  Problem  start  time  if  other  than  zero  may  be  specified  as 
TIMEO.  Other  control  constants  used  or  activated  are  TIMEN,  TIMEM,  CSGMIN, 
CSGFAC,  DTIMEU,  DTIMEL,  DTIMEH,  DTMPCA,  DTMPCC,  ATMPCA,  ATMPCC, 
ARLXCA,  ARLXCC,  DAMPA,  OPEITR,  BACKUP,  LINECT,  and  PAGECT. 


Calling  Sequence— CNEXPN  —  This  subroutine  utilizes  one  dynamic  storage  core 
location  for  each  diffusion  and  arithmetic  node. 
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Execution  Subroutine  CNFWBK 

Purpose— This  subroutine  performs  transient  thermal  analysis  by  implicit  forward- 
backward  differencing.  The  LPCS  option  is  required  and  allows  the  simultaneous  set  of 
equations  to  be  solved  by  “successive  point”  iterations.  During  the  first  iteration  for  a 
time  step,  the  capacitance  values  are  doubled  and  divided  by  the  time  step  and  the  energy 
transfer  rates  based  on  old  temperatures  are  added  to  the  source  locations.  Upon  com¬ 
pleting  the  time  step  the  capacitance  values  are  returned  to  their  original  state.  The  itera¬ 
tion  looping,  convergence  criteria  and  other  control  constant  checks  are  identical  to 
CNBACK.  The  time  step  checks  and  calculations  and  look  ahead  feature  are  identical  to 
that  used  for  CNBACK. 

The  automatic  radiation  transfer  damping  and  extrapolation  method  of  acceleration 
mentioned  under  the  r'NBACK  subroutine  writeup  are  also  employed  in  this  subroutine. 
Diffusion  and/or  arithmetic  temperature  calculations  may  be  damped  through  use  of 
DAMPD  and/or  DAMPA  respectively.  Control  constants  BACKUP  and  OPEITR  are  con¬ 
tinuously  checked.  CNFWBK  internally  performs  forward-backward  differencing  of 
boundary  conditions.  For  this  reason  the  user  should  utilize  TIMEN  as  the  appropriate 
independent  variable  in  Variables  1  operations. 

It  is  interesting  to  note  that  CNFWBK  generally  converges  in  25%  fewer  iterations 
than  CNBACK.  The  probable  reason  for  this  is  that  the  boundary  of  the  mathematical 
system  is  better  defined.  While  every  future  temperature  node  under  CNBACK  is  con¬ 
nected  to  its  present  temperature,  under  CNFWBK  every  future  temperature  node  is  also 
receiving  an  impressed  source  based  on  the  present  temperature. 


Restrictions— The  LPCS  option  is  required.  Control  constants  TIMEND,  OUTPUT, 
DTIMEI,  NLOOP  and  DRLXCA  and/or  ARLXCA  must  be  specified.  Other  control  con¬ 
stants  which  are  used  or  activated  are  TIMEN,  TIMEO,  TIMEM,  CSGMIN,  DTIMEU, 
DTIMEH,  DTMPCA,  DTMPCC,  ATMPCA,  ATMPCC,  DAMPD,  DAMPA,  DRLXCC  and/or 
ARLXCC,  LOOPCT,  BACKUP,  OPEITR,  LINECT,  and  PAGECT. 


Calling  Sequence  -CNFWBK  -  Tills  subroutine  utilizes  three  dynamic  storage  core 
locations  for  each  diffusion  node  and  one  for  each  arithmetic  and  boundary  node. 


Execution  Subroutine  CNBACK 

Purpose— This  subroutine  performs  transient  thermal  analysis  by  implicit  backward 
differencing.  The  LPCS  option  is  required  and  allows  the  simultaneous  set  of  equations 
to  be  solved  by  “successive  point”  iteration.  Each  third  iteration,  diffusion  node  tempera¬ 
tures  which  trace  a  continuous  decreasing  slope  receive  an  extrapolation  on  their  error 
function  curve  in  an  attempt  to  accelerate  convergence.  For  convergence  criteria  the  user 
is  required  to  specify  NLOOP  and  DRLXCA  and/or  ARLXCA.  If  the  number  of  itera¬ 
tions  during  a  time  step  exceeds  NLOOP  a  message  is  printed  but  the  problem  proceeds. 
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Variables  1  is  performed  only  once  for  each  time  step.  Since  this  subroutine  is  im¬ 
plicit  the  user  must  specify  the  time  step  to  be  us' d  as  DTIMEI  in  addition  to  TIMEND 
and  OUTPUT.  The  look  ahead  feature  for  the  time  step  calculation  in  CNFRWD  is  used 
as  are  the  checks  for  DTIMEH,  DTMPCA  and  ATMPCA  but  not  DTIMEL.  Damping  of 
the  solutions  can  be  achieved  through  use  of  control  constants  DAMPD  and/or  DAMPA. 
Control  constants  BACKUP  and  OPEITR  are  continuously  checked. 

Implicit  methods  of  solution  often  oscillate  at  start  up  or  for  boundary  step  changes 
when  radiation  conductors  are  present.  CNBACK  contains  an  automatic  damping  feature 
which  is  applied  to  radiation  conductors.  The  radiation  transfer  to  a  node  is  calculated 
for  its  present  temperature  and  a  temporary  new  temperature  is  calculated.  Then  the 
radiation  transfer  is  recalculated  and  the  final  node  temperature  is  calculated  based  on 
the  arithmetic  mean  of  the  two  radiation  transfer  calculations.  This  automatic  radiation 
damping  has  proven  to  be  quite  successful  and  lessens  the  need  for  use  of  DAMPD  and 
DAMPA. 


Restrictions— The  LPCS  option  is  required.  Control  constants  TIMEND,  OUTPUT, 
DTIMEI,  NLOOP  and  DRLXCA  and/or  ARLXCA  must  be  specified.  Other  control  con¬ 
stants  which  are  used  on  activated  are:  TIMEN,  TIMEO,  TIMEM,  CSGMIN,  DTIMEV, 
DTIMEH,  DTMPCA,  DTMPCC,  ATMPCA,  ATMPCC,  DAMPD,  DAMPA,  DRLXCC  and/or 
ARLXCC,  LOOPCT,  BACKUP,  OPEITR,  LINECT,  and  PAGECT. 


Calling  Sequence— CNBACK  —  This  subroutine  utilizes  three  dynamic  storage  core 
locations  for  each  diffusion  node  and  one  for  each  arithmetic  and  boundary  node. 
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Interpolation  Subroutines 


Name  Page 
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D11DIM,  D11MDI .  60 
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CVQ1HT,  CVQ1WM,  GSL0PE,  PSINTR,  PSNTWM .  62 
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Subroutine  LAGRAN  or  LGRNDA 

Purpose— These  subroutines  perform  Lagrangian  interpolation  of  up  to  order  50.  The 
first  requires  one  doublet  array  of  X,  Y  pairs  while  the  second  requires  two  singlet  arrays, 
one  of  X’s  and  the  other  of  Y’s.  They  contain  an  extrapolation  feature  such  that  if  the 
X  value  falls  outside  the  range  of  the  independent  variable  the  nearest  dependent  Y  vari¬ 
able  value  is  returned  and  no  error  is  noted. 

n  ^  X  X* 

Y  =  Pn<X>  =  Z  Yk  II  x~Tx-  *  r*  =  1.2, 3,-.., 50 max. 
k  =  0  i=0  k  ' 

i*k 

Restrictions— All  values  must  be  floating  point  except  N  which  is  the  ord-r  of  inter¬ 
polation  plus  one  and  must  be  an  integer.  The  independent  variable  values  must  be  in 
ascending  order. 

Calling  Sequence- LAGRAN(X,Y,A(IC),N)  or  LGRNDA(X,Y,AX(1C),AY(IC),N) 

NOTE:  A  doublet  array  is  formed  as  follows: 

IC,  XI,  Yl,  X2,  Y2,  X3,  Y3,...,  XN,  YN 
where  IC  =  2*N  (set  by  program). 

Singlet  arrays  are  formed  as  follows: 

IC,  XI,  X2,  X3,  ....  XN 
IC,  Yl,  Y2,  Y3, ...,  YN 
and  IC  =  N  (set  by  program). 

Subroutine  D1DEG1  or  D1D1DA 

Purpose— These  subroutines  perform  single  variable  linear  interpolation  on  doublet  or 
singlet  arrays  respectively.  They  are  self-contained  subroutines  that  are  called  upon  by 
virtually  all  other  linear  interpolation  subroutines. 

Restrictions— All  values  must  be  floating  point  numbers.  The  X  independent  variable 
values  must  be  in  ascending  order. 

Calling  Sequence- D1DEG1(X,A(IC),Y)  or  D1D1DA(X,AX(1C),AY(IC),Y) 

Subroutine  D1D1WM  or  D11MDA 

Purpose— These  subroutines  perform  single-variable,  linear  interpolation  by  calling  on 
D1DEG1  or  D1D1DA,  respectively.  However,  the  interpolated  answer  is  multiplied  by 
the  value  addressed  as  Z  prior  to  being  returned  as  Y. 
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Restrictions— Same  as  D1DEG1  or  D1D1DA,  and  Z  must  be  a  floating-point  number. 
Calling  Sequence— D1D1WM(X,A(1C),Z,Y)  or  D11MDA(X>AX(IC),AY(IC),Z,Y) 


Subroutine  D1MDG1  or  D1M1DA 

Purpose— These  subroutines  use  the  arithmetic  mean  of  two  input  values  as  the  inde¬ 
pendent  variable  for  linear  interpolation.  They  require  a  doublet  or  two  singlet  arrays, 
respectively. 

Restrictions— See  D1DEG1  or  D1D1DA  as  they  are  called  on,  respectively. 

Calling  Sequence— D1MDGJ  1X1,X2,A(1C),Y)  or 

D1M1DA(X1,X2,AX(IC),AY(IC),Y) 


Subroutine  D1M1WM  or  D1M1MD 

Purpose— These  subroutines  use  the  arithmetic  mean  of  two  input  values  as  the  inde¬ 
pendent  variable  for  linear  interpolation.  The  interpolated  answer  is  multiplied  by  the  Z 
value  prior  to  being  returned  as  Y. 

Restrictions— Same  as  D1MDG1  or  D1M1DA,  and  Z  must  be  a  floating-point  number. 

Calling  Sequence— D1M1WM(X1,X2,A(IC),Z,Y)  or 

D1M1MD(X1,X2,AX(IC),AY(!C),Z,Y) 


Subroutine  D1DEG2  or  D1D2DA 

Purpose— These  subroutines  perform  single-variable  parabolic  interpolation.  The  first 
requires  a  doublet  array  of  X,  Y  pairs  while  the  second  requires  singlet  arrays  of  X  and 
Y  values.  They  call  on  subroutines  LAGRAN  and  LGRNDA,  respectively. 

Restrictions— See  LAGRAN  or  LGRNDA. 

Calling  Sequence -D1DEG2(X,A(IC),Y)  or  D1D2DA(X,AX(IC),AY(IC),Y) 


Subroutine  D1D2WN1  or  D12MDA 

Purpose— These  subroutines  perform  single-variable  parabolic  interpolation  by  calling 
on  LAGRAN  or  LGRNDA,  respectively.  However,  the  interpolated  answer  is  multiplied 
by  the  value  addressed  as  Z  prior  to  being  returned  as  Y. 

Restrictions— Same  as  LAGRAN  or  LGRNDA,  and  Z  must  be  a  floating-point 
number. 
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Calling  Sequence— D1D2WM(X,A(IC),Z,Y)  or 

D12MDA(X,AX(IC),AY(IC),Z,Y) 


Subroutine  D1MDG2  or  D1M2DA 

Purpose— These  subroutines  use  the  arithmetic  mean  of  two  input  values  as  the  inde¬ 
pendent  variable  for  parabolic  interpolation.  They  require  a  doublet  or  two  singlet  arrays, 
respectively. 

Restrictions— See  LAGRAN  or  LGRNDA  as  they  are  called. 

Calling  Sequence— D1MDG2(X1,X2,A(IC),Y)  or 

D1M2DA(X1,X2,AX(IC),AY(IC),Y) 


Subroutine  D1M2WM  or  D1M2MD 

Purpose— These  subroutines  use  the  arithmetic  mean  of  two  input  values  as  the  inde¬ 
pendent  variable  for  parabolic  interpolation.  The  interpolated  answer  is  multiplied  by  the 
Z  value  prior  to  being  returned  as  Y. 

Restrictions— Same  as  D1MDG2  or  D1M2DA,  and  Z  must  be  a  floating  point  number. 

Calling  Sequence— D1M2WM(X1,X2,A{IC),Z,Y)  or 

D1M2MD(X1,X2,AX(IC),AY(IC),Z,Y) 


Subroutine  D1DG1I  or  D1D1IM  or  DIDIMi 

Purpose— These  subroutines  perform  single-variable  linear  interpolation  on  an  array 
of  X’s  to  obtain  an  array  of  Y’s.  D1D1IM  multiplies  all  interpolated  values  by  a  constant 
Z  value  while  DIDIMI  allows  a  unique  Z  value  for  each  X  value.  They  all  call  on 
D1DEG1. 

Restrictions— The  number  of  input  X’s  must  be  supplied  as  the  integer  N  and  agree 
with  the  number  of  Y  and  Z  locations  where  applicable.  Z  values  must  be  floating-point 
numbers. 

Calling  Sequence— DIDG1I(N,X(DV),A(IC),Y(DV))  or 
D1D1IM(N,X(DV),A(IC),Z,Y(DV))  or 
D1D1MI(N,X(DV),A(IC),Z(DV),Y(DV)) 


Subroutine  D11DAI  or  D11DIM  or  D11MDI 

Purpose— These  subroutines  are  virtually  identical  to  D1DG1I,  D1D1IM,  and  DIDIMI, 
respectively.  The  difference  is  that  they  require  singlet  arrays  for  interpolation  and  call 
on  D1D1DA. 
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Restrictions— Same  as  D1DG1I,  D1D1IM,  and  D1D1MI. 

Calling  Sequence- D11DAI(N,X(DV),AX(1C),AY(IC),Y(DV))  or 
D11D1M(N,X(DV),AX(IC),AY(IC),Z,Y(DV))  or 
D11MDI(N,X(DV),AX(IC),AY(IC),Z(DV),Y(DV)) 


Subroutine  D11CYL  or  DAI  ICY 

Purpose— These  subroutines  reduce  core  storage  requirements  for  cyclical  interpola¬ 
tion  arrays.  The  arrays  need  cover  one  period  only,  and  the  period  (PR)  must  be  speci¬ 
fied  as  the  first  argument.  Linear  interpolation  is  performed,  and  the  independent  variable 
must  be  in  ascending  order. 

Restrictions— All  values  must  be  floating  point.  Subroutine  INTRFC  is  called  on  by 
both  D11CYL  and  DA11CY,  then  D1DEG1  or  D1D1DA,  respectively. 

Calling  Sequence- D11CYL(PR,X,A(IC),Y)  or 

DA11CY(PR,X,AX(IC),AY(IC),Y) 


Subroutine  D12CYL  or  DA12CY 

Purpose— These  subroutines  are  virtually  identical  to  D11CYL  and  DAI  ICY,  except 
that  parabolic  interpolation  is  performed. 

Restrictions— See  D11CYL  and  DA11CY.  Subroutines  LAGRAN  and  LGRNDA, 
respectively,  are  called  on. 

Calling  Sequence- D12CYL(PR,X,A(IC),Y)  or  DA12CY(PR,X,AX(1C),AY(IC),Y) 


Subroutine  D11MCY  or  DA11MC 

Purpose— These  subroutines  are  virtually  identical  to  D11CYL  or  DA11CY,  except 
that  the  interpolated  answer  is  multiplied  by  the  floating-point  Z  value  prior  to  being 
returned  as  Y. 

Restrictions— Call  on  subroutines  D1DEG1  and  D1D1DA,  respectively. 

Calling  Sequence— D11MCY(PR,X,A(IC),Z,Y)  or 

DA11MC(PR,X,AX(IC),AY(IC),Z,Y) 


Subroutine  D12MCY  or  DA12MC 

Purpose- -These  subroutines  are  virtually  identical  to  D11MCY  and  DA11MC  except 
that  parabolic  interpolation  is  performed. 
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Restrictions— C alls  on  subroutines  LAGRAN  and  LGRNDA,  respectively. 

Calling  Sequence— D12MCY(PR,X,A(IC),Z,Y)  or 

DA12MC(PR,X,AX(IC),AY(IC),Z,Y) 


Subroutine  CVQ1HT  or  CVQ1WM 

Purpose— These  subroutines  perform  two  single-variable  linear  interpolations.  The 
interpolation  arrays  must  have  the  same  independent  variable  X  and  dependent  variables 
of,  say,  R(X)  and  S(X).  Additional  arguments  of  Y,  Z,  and  T  complete  the  data  values. 
The  postinterpolation  calculations  are,  respectively: 

Y  =  S(X)*(R(X)-T)  or 

Y  =  Z*S(X){R(X)-T). 

Restrictions— Interpolation  arrays  must  be  of  the  doublet  type  and  have  a  common 
independent  variable.  All  values  must  be  floating-point  numbers. 

Calling  Sequence— CVQ1HT(X,AR(IC),AS(IC),T,Y)  or 

CVQ1WM(X,AR(IC),AS(IC),T,Z,Y) 


Subroutine  GSLOPE 

Purpose— This  subroutine  will  generate  a  slope  array  so  that  point  slope  interpolation 
subroutines  can  be  used  instead  of  standard  linear  interpolation  subroutines.  The  user 
must  address  two  singlet  arrays,  and  a  singlet  slope  array  will  be  produced. 

Restrictions— The  X  independent-variable  array  must  be  in  ascending  order.  All  ar¬ 
rays  must  be  of  equal  length  and  contain  floating-point  numbers. 

Calling  Sequence-G SLOPE(AX(IC),AY(IC),AS(IC)) 


Subroutine  PSINTR  or  PSNTWM 

Purpose— These  subroutines  perform  linear  interpolation  and  require  arrays  of  the  Y 
points  and  slopes  which  correspond  to  the  independent  variable  X  array.  All  values  must 
be  floating-point  numbers.  PSNTWM  multiplies  the  interpolated  answer  by  Z  prior  to 
returning  it  as  Y. 

Restrictions— The  independent  X  and  dependent  Y  and  slope  arrays  must  be  of 
equal  length. 

Calling  Sequence- PSINTR(X,AX(IC),AY(IC),AS(IC),Y)  or 
PSNTWM(X,AX(IC),AY(IC),AS(IC),Z,Y) 
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Bivariate  Array  Format,  Z  =  f(X,Y) 

Bivariate  arrays  must  be  rectangular  and  full  and  must  be  entered  in  the  following 
row  order: 


IC,N  ,X  1,X  2,X  3,  . . . ,  X  N 
Yl, Zll, Z12, Z13,  ....  Z1N 
Y2,  Z21,  Z22,  Z23,  . .  . ,  Z2N 


YM,ZM1,ZM2,ZM3,  . . .  ,ZMN 

where  N  is  the  integer  number  of  X  variables.  All  other  values  must  be  floating-point  num¬ 
bers,  and  the  X  and  Y  values  must  be  in  ascending  order. 


Subroutine  BVSPSA  or  BVSPDA 

Purpose— These  subroutines  use  an  input  Y  argument  to  address  a  bivariate  array  and 
pull  off  a  singlet  array  of  Z’s  correspond:ng  to  the  X’s  or  pull  off  a  doublet  array  of  X,Z 
values,  respectively.  The  integer  count  for  the  constructed  arrays  must  be  exactly  N  or  2*N, 
respectively.  To  use  the  singlet  array  for  an  interpolation  call,  reach  the  X  array  by  address¬ 
ing  the  N  in  the  bivariate  array. 

Restrictions— As  stated  above,  and  all  values  must  be  floating  point. 

Calling  Sequence— BVSPSA(Y,BA(IC),AZ(JC))  or  BVSPDA(Y,BA(IC),AXZ(IC)) 


Subroutine  D2DEG1  or  D2DEG2 

Purpose— These  subroutines  perform  bivariate  linear  and  parabolic  interpolation,  re¬ 
spectively.  The  arrays  must  be  formatted  as  shown  for  Bivariate  Array  Format. 

Restrictions— For  D2DEG1,  N  >  2,  M  ^  2l  See  bivariate 
For  D2DEG2,  N  ^  3,  M  >  3j  array  format 

Cu//mgSeq::ence-D2DEGl(X,Y,DA(lC),Z)  or  D2DEG2(X,Y,BA(1C),Z) 


Subroutine  D2D1WM  or  D2D2WM 

Purpose— These  subroutines  perform  bivariate  linear  or  parabolic  interpolation  by 
calling  on  D2DEG1  or  D2DEG2,  respectively.  The  interpolated  answer  is  multiplied  by 
the  W  value  prior  to  being  returned  as  Z. 

Restrictions— Hame  as  D2DEG1  or  D2DEG2,  and  W  must  be  a  floating-point  value. 
Calling  Sequence — D2D 1  WM(X. Y,BA(  1C),W,Z)  or  D2D2WM(X.Y,BA{IC),\V,Z) 
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MAn  I  t,  ULAL  t  Gy,  £  '  •  .£**>  „  •  \ 

Subroutine  D2MXD1  or  D2MXD2  '  s>>  <■'■'" i'  ■  '■  s* 

Purpose— These  subroutines  are  virtually  identical  to  D2D3SGA  -  and  D ?. D  '  ■ 

that  the  arithmetic  mean  of  two  X  values  is  used  as  the',  A -independent  variable  for. .  it  ,  .Y,.-’  .j," 
interpolation.  \  t  ..  '  j._, 

'3>,  '*  e  '  ...  .  I  I 

Restrictions— Same  as  D2DEG1  or  D2DKG2.  -.  V  " 


Y 


Calling  Scgue/ice— D2MXD1  (XI, X2,Y,BA(1C),Z)  or  D2M XD2(X  1  ,X2,  Y.-ip  A I IC }.Z JU-; 


Subroutine  D2MX1M  or  D2MX2M  . < 

v  ■.* 

Purpose— These  subroutines  are  virtually  identical  to  U2D1WM  and  02D2WM j:xcepi 
that  the  arithmetic  mean  of  two  X  values  is  used  as  the  X-iridependenl  variable  .for  j 
interpolation.  ,  N  .1;  ,  ‘ 

‘  ,  ft”." 

Restrictions— Same  as,  D2D1WM  and  D202WM.  "  "$A*. 

'  '• 

Calling  Sequence— D2MX1M(X1.X2,Y,BA(1C),W,Z)  or  -  •  "A  * 

D2MX2M(X1,X2,Y,BA(IC),W,Z) 


.  ’If'*'. 


Trivariate  Array  Format,  T  =  f(X,Y,Z)  -> 

Trivariate  anrays  may  be  thought  of  as  two  or  more  bivariate  arrays,  .each  to  vans  a* 
array  a  function  of  a  third  independent  variab’e  Z.  Trivariate  arrays  must,  be  entered  >r,v 
row  order  and  be  constructed  as  follows:  *  <  t. 

-  4'. 

IC,NX1,NY1,  Z1,X  1,X  2,X  3, _ X  N 


NX1.NY1,  Z1,X  1,X  2,  X  3,  . 

.  .  .  X  N 

Yl,  Til,  T12,  T13,  . 

. . ,  TIN 

Y2,  T21,  T22,  T23,  . 

.  .  ,  T2N 

YM,TM1,TM2,TM3,  . 

. . ,TMN 

NX2,NY2,  Z2,  X  1,  X  2,  X  3,  . 

...  X  J 

Yl,  Til,  T1 2, T13,  . 

.  .  ,  T1J 

Y2,  T21.T22,  T23,  . 

.  .  ,  T2J 

YK,TK1,TK2,TK3,  . 

.  .  ,  TKJ 

NX3.NY3,  Z3,  .  .  . 

...  . 

...  ‘  V  *•; 

1  .?/  •  ■ 

' 1  /  J/, 

i  * 

7  •  : 


v  %v  y  v; 

*  •  i  >  /- ' 


The  trivariate  array  may  consist  of  as  many  bivariate  “sheets”  as  desired  The  dum¬ 
ber  of  X  and  Y  values  in  each  sheet  must  be  specified  as  integers  (NX -NY ).  The  ’Sheets" 
must  be  rectangular  and  full  but  need  not  be  identical  in  size. 
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Subroutine  D3DEG1  or  D3D1WM 

Purpose— These  subroutines  perform  trivariate  linear  interpolation.  The  interpolation 
array  must  be  constructed  as  shown  for  the  Trivariate  Array  Format.  Subroutine  D2DEG1 
is  called  on,  which  calls  on  D1DEG1.  Hence,  the  linear  extrapolation  feature  of  these 
routines  applies.  Subroutine  D3D1WM  multiplies  the  interpolated  answer  by  F  prior  to 
returning  it  as  T. 

Restrictions— See  Trivariate  Array  Format.  F  must  be  a  floating-point  value. 

Calling  Sequence— D3DEGl(X,Y.Z,TA(IC),T)  or  D3D1WM(X,Y,Z,TA(IC),F,T) 


Subroutine  VARCSM  or  VARCCM  or  VARC1  or  VARC2 

Purpose— These  are  linear  interpolation  subroutines  which  are  set  up  as  Variables  1 
calls  by  the  preprocessor  when  processing  the  CGS  and  CGD  mnemonic  codes  in  the  nodal 
data  block.  VARCSM  is  utilized  for  the  CGS  code.  VARCCM  uinu’eu  for  the  CGD 
code  when  two  array  arguments  appear.  VARC1  and  VARC2  are  used  for  the  CGD  code 
when  either  the  first  or  second  respective  array  arguments  are  entered  as  a  constant.  The 
following  mnemonic  codes  in  the  nodal  block 

8 

I 

4 

CGS  I ,80. , A1 , 10.2 
CGD  2,80. , A 1 , 10.2,A2, 1 .6 
CGD  3, 80., I. 4, 5. I, A2, 1.6 
CGD  4, 80. , At, b. I, 6. 3, 3. 7 

would  cause  the  construction  in  Variables  1  of 

12 

1 

VARCSMtTI ,C1 , A) ,10.2) 

VARCCM ( f2,C2,AI , 1 0.2,  A2, 1 .6) 

VARCI (T3,C3,l.4,b.l ,A2, 1.6) 

VARC2(T4,C4,AI ,S. 1 ,6.  3, 8,  7) 

The  second  call  causes  the  sum  of  two  interpolations  with  multiplications  to  be  used  as 
the  C2  value.  The  latter  two  calls  only  perform  one  interpolation  but  use  the  sum  of  the 
two  products  as  the  C  value. 

Restrictions— The  array  arguments  must  address  the  integer  count. 

Calling  Sequence — V A RCSM(T.C, At IC ),F )  or  VARCCM(T,C,A1(IC),F1,A2(1C),F2)  or 
VARCI (T,C,X,FI,A2( 1C). F2)  or  VARC2(T,C,A1(IC),F1  .X.F2) 
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Subroutine  VARGSM  or  VARGCM  or  VARG1  or  VARG2 

Purpose— These  are  linear  interpolation  subroutines  set  up  as  Variables  1  calls  by  the 
preprocessor  when  processing  the  CGS  and  CGD  mnemonic  codes  in  the  conductor  data 
block.  They  are  similar  to  the  preceding  four  calls  for  the  nodal  data  block  except  that 
the  conductor  argument  is  first  followed  by  two  temperature  arguments.  VARGSM  is 
used  for  the  CGS  code.  If  the  F  value  is  positive,  the  mean  of  the  two  addressed  tem¬ 
peratures  is  used  for  interpolation.  If  it  is  negative,  only  T1  is  used  for  interpolation  and 
the  absolute  value  of  F  is  used  as  a  multiplier.  VAROCM,  VARG1,  and  VARG2  per¬ 
form  the  one  or  two  interpolations  required,  multiply  by  the  F  values  to  obtain  G1  and 
G2  components,  and  then  calculate  G  as 

G  =  1.0/(1.0/G1  +  1.0/G2). 

Restrictions— The  array  arguments  must  address  the  integer  count. 

Calling  Sequence— VARGSM(G,T1,T2,A(IC),F)  or 

VARGCM(G,T1,T2,A1(1C),F1,A2(1C),F2)  or 
VARG1(G.T1,T2,X,F1,A2(IC),F2)  or 
VARG2(G,T1,T2,A1(IC),F1,X,F2) 
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Arithmetic  Subroutines 


Name  Page 

FL0AT,  FIX,  INTRFC,  SHFTV,  SHFTVR,  FLIP, 

SETPLS,  ARYPLS  .  68 

SETMNS,  ARYMNS,  ADD,  ADDFIX,  ADDARY,  ARYADD  .  69 

SUB,  SUBFIX,  SUBARY,  ARYSUB,  MLTPLY,  MPYFIX, 

MPYARY,  ARYMPY .  70 

DIVIDE,  DIVFIX,  DIVARY,  ARYDIV,  GENARY  .  71 

BLDARY,  BRKARY,  BKARAD,  STFSEP,  SCALE,  .  72 

STFSEQ,  STFSQS,  SUMARY,  MAXDAR,  MXDRAL  .  73 

ARYINV,  ARINDV,  ADDINV.  ADARIN,  ST0ARY,  ARYST0  .  74 

SCLDEP,  SCLIND,  SLDARY,  SLDARD,  SPLIT,  J0IN .  75 

SPREAD,  QMETER,  RDTNQS,  QMTRI,  QF0RCE,  QINTEG, 

QINTGI  .  76 

CINSIN,  BINARY,  CINC0S,  C0SARY,  CINTAN,  TANARY, 

ARCSIN,  ASNARY  .  77 

ARC0S,  ACSARY,  ARCTAN,  ATNARY,  EXPNTL, 

ARYEXP,  EXPARY  .  78 

L0GT,  L0GTAR,  L0GE,  L0GEAR,  SQR00T,  SQR0TI,  CMPXSR, 

CSQRI  .  79 

CMPXMP,  CMPYI,  CMPXDV,  CDIVI,  NEWTRT,  NEWRT4 .  80 

PLYNML,  PLYARY,  SMPINT,  TRPZD,  TRPZDA  .  81 

PRESS,  SPR.ESS,  EFFG,  EFFEMS  .  82 
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Subroutine  FLOAT  or  FIX  or  INTRFC 

Purpose— Subroutine  FLOAT  will  convert  an  integer  to  a  floating-point  number. 
Subroutine  FIX  will  convert  a  floating-point  number  to  an  integer.  Subroutine  INTRFC 
will  fracture  a  floating-point  number  to  yield  the  largest  integer  value  possible  and  the 
remainder  or  fractional  portion  is  a  floating-point  number.  Their  respective  operations 
are 


X  =  N 
orN  =  X 
or  N  =  X 
Y  =  N 
F  =  X-Y 

Restrictions— X  and  F  arguments  must  address  floating-point  values  and  the  N  argu¬ 
ment  must  address  an  integer. 

Calling  Sequence- FLOAT(N,X)  or  FIX(X,N)  or  INTRFC(X,N,F) 


Subroutine  SIIFTV  or  SHFTVR  or  FLIP 

Purpose— Subroutine  SHFTV  will  shift  a  sequence  of  data  from  one  array  to  another. 
Subroutine  SHFTVR  will  shift  a  sequence  of  data  from  one  array  and  place  it  in  another 
array  in  reverse  order.  Subroutine  FLIP  will  reverse  an  array  in  its  own  array  location. 
Their  respective  operations  are 

A(i>  “  B(i),  i  =  1,  N 

or  A(N-i+l)  =  B(i),  i  =  1,  N 

°r  A(i,new  =  A(N-i+2)0y,  i  =  2,  N+l. 

The  answer  array  may  not  be  overlayed  into  the  input  array. 

Restrictions— The  data  values  to  be  shifted  or  reversed  in  order  may  be  anything. 

The  N  must  be  an  integer. 

Calling  Sequence— SHFTV(N,B(DV),A(DV))  or  SHFTVR(N,B(DV),A(DV))  or 
FLIP(AOO) 


Subroutine  SETPLS  or  ARYPLS 

Purpose— SETPLS  will  set  the  sign  positive  for  a  variable  number  of  arguments  while 
ARYPLS  will  set  the  sign  positive  for  every  data  value  in  a  specified  length  array. 

Restrictions— The  values  addressed  may  be  either  integers  or  floating-point  numbers. 
The  number  (N)  of  data  values  in  the  array  must  be  specified  as  an  integer. 
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Calling  Sequence — SETPLS(  A,  B,C...)  or  ARYPLS(N,A(DV)) 

where  N  may  be  a  literial  integer  or  the  address  of  a  location  containing  an  integer,  and 
A(DV)  addresses  the  first  data  value  in  the  array. 


Subroutine  SETMNS  or  ARYMNS 

Purpose— SETMNS  will  set  the  sign  negative  for  a  variable  number  of  arguments, 
while  ARYMNS  will  set  the  sign  negative  for  every  data  value  in  a  specified  length  array. 

Restrictions-  '’’he  values  addressed  may  be  either  integers  or  floating-point  numbers. 
The  number  (N)  of  data  values  in  the  array  must  be  specified  as  an  integer. 

Calling  Sequence— SETMNS(A,B,C, . . . )  or  ARYMNS(N,A(DV)) 

where  N  may  be  a  literial  integer  or  the  address  of  a  location  containing  an  integer  and 
A(DV)  addresses  the  first  data  value  in  the  array. 


Subroutine  ADD  or  ADDFIX 

Purpose— To  sum  a  variable  number  of  floating-point  or  integer  numbers,  respectively. 
S  =  SXj,  i  =  1,  2,  3,...,  N,  N  >  2 

Restrictions— Subroutine  ADD  is  for  floating-point  numbers,  while  subroutine 
ADDFIX  is  for  integers. 

Calling  Sequence- ADD(X1,X2,X3, . . .  ,XN,S)  or  ADDFIX(X1,X2,X3, . . .  ,XN,S) 


Subroutine  ADDARY  or  ARY  ADD 

Purpose— Subroutine  ADDARY  will  add  the  corresponding  elements  of  two  specified 
length  arrays  to  form  a  third  array.  Subroutine  ARY  ADD  will  add  a  constant  value  to 
every  element  in  an  array  to  form  a  new  array.  Their  respective  operations  are 

Aj  =  B;  +  C„  i  =  1,  N 
or  A;  =  Bj  +  C,  i  =  1,  N. 

The  answer  array  may  be  overlayed  into  one  of  the  input  array  areas. 

Restrictions— All  data  values  to  be  operated  on  must  be  floating-point  numbers.  The 
array  length  N  must  be  an  integer. 

Calling  Sequence— ADDARY(N,B(DV),C(DV),A(DV))  or 
ARYADD(N,B(DV),C,A(DV)) 
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Subroutine  SUB  or  SUBFIX 

Purpose— To  subtract  a  variable  number  of  floating-point  or  integer  numoers, 
respectively, 

R  =  Y  -  IX;,  i  =  1,  2,  3,  N,  N  >  1 

Restrictions— Subroutine  SUB  is  for  floating-point  numbers  while  subroutine  SUBFIX 
is  for  integers. 

Calling  Sequence -SUB(  Y,X1,X2,X3,. . .  ,XN,R)  or 
SUBF1X(Y,X1,X2,X3 . XN,R) 


Subroutine  SUBARY  or  ARYSUB 

Purpose— Subroutine  SUBARY  will  subtract  the  corresponding  elements  of  one  array 
from  another  to  form  a  third  array.  Subroutine  ARYSUB  will  subtract  a  constant  value 
from  every  element  in  an  array  to  form  a  new  array.  Their  respective  operations  are 

Aj  =  Bj  -  Cj,  i  -  1,  N 
or  Aj  =  Bj  -  C,  i  =  1,  N 

The  answer  array  may  be  overlayed  into  one  of  the  input  array  areas. 

Restrictions— All  data  values  to  be  operated  on  must  be  floating-point  numbers.  The 
array  length  N  must  be  an  integer. 

Calling  Sequence --SUBARY(N,B(DV),C(DV),A(DV))  or 
ARYSUB(N,B(DV),C,A(DV)) 


Subroutine  MLTPLY  or  MPYFIX 

Purpose— To  multiply  a  variable  number  of  floating-point  or  integer  numbers, 
respectively. 

P  »  XI  *  X2  *  X3  *  ...  *  XN,  U>2 

Restrictions— Subroutine  MLTPLY  is  for  floating-point  numbers,  while  subroutine 
MPYFIX  is  for  integers. 

Calling  Sequence — M LTPLY ( X 1  ,X2,X3 . XN,P)  or  MPYFIX(X1,X2,X3 . XN,P) 


Subroutine  MPYARY  or  ARYMPY 

Purpose— Subroutine  MPYARY  will  multiply  the  corresponding  elements  of  two 
arrays  to  form  a  third.  Subroutine  ARYMPY  will  multiply  a  constant  value  times  each 
element  of  an  array  to  form  a  new  array.  Their  respective  operations  are 
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A;  =  Bj  *  Cj,  i  =  1,  N 
or  Aj  =  Bj  *  C,  i  =  1,  N 

The  answer  array  may  be  overlayed  into  one  of  the  input  array  areas. 

Restrictions— All  data  values  to  be  operated  on  must  be  floating-point  numbers.  The 
array  length  N  must  be  an  integer. 

Calling  Sequence- MPYARY(N,B(DV),C(DV),A(DV))  or 

ARYMPY(N,B(DV),C,A(DV)) 


Subroutine  DIVIDE  or  DIVFIX 

Purpose— These  subroutines  are  used  to  perform  a  division  of  floating-point  or  inte¬ 
ger  numbers,  respectively; 

Q  =  Y/£  Xj,  i  =  1,  2,3 . N,  N>1. 

Restrictions— Subroutine  DIVIDE  is  for  floating-point  numbers,  while  DIVFIX  is  for 
integers. 

Calling  Sequence- DIVIDE(Y,X1,X2,X3, . . .  ,XN,Q)  or 
DIVFIX(Y,X1,X2,X3, . . .  ,XN,Q) 


Subroutine  DIVARY  or  ARYDIV 

Purpose— Subroutine  DIVARY  will  divide  the  elements  of  one  array  into  the  corre¬ 
sponding  elements  of  another  array  to  produce  a  third  array.  Subroutine  ARYDIV  will 
divide  each  element  of  an  array  by  a  constant  value  to  produce  a  new  array.  Their  re¬ 
spective  operations  are 

Aj  »  Bj/Cj,  i  =  1,  N 
or  Aj  =  Bi/C,  i  =  1,  N. 

The  answer  array  may  be  overlayed  into  one  of  the  input  array  areas. 

Restrictions— All  data  values  to  be  opei-I^d  on  must  be  floating-point  numbers.  The 
array  length  N  must  be  an  integer. 

Calling  Sequence— D1VARY(N,B(DV),C(DV),A(DV))  or 
ARYDIV(N,B(DV),C,A(DV)) 


Subroutine  GENARY 

Purpose— This  subroutine  will  generate  an  array  of  equally  incremented  ascending 
values.  The  user  must  supply  the  minimum  value,  maximum  value,  number  of  values  in 
the  array  to  be  generated,  and  the  space  for  the  generated  array. 
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Restrictions— All  numbers  must  be  floating  point. 
Calling  Seque«cc~GENARY(B(DV),A(DV)) 


whe»e 


B(l)  =  minimum  value 
B(2)  =  maximum  value 

B(3)  =  length  of  array  to  be  generated  (floating  point). 


Subroutine  BLDARY 

Purpose— This  subroutine  will  build  an  array  from  a  variable  number  of  arguments  in 
the  order  listed.  The  operation  performed  is 

A;  =  Xj,  i  =  1,  N. 

Restrictions— Data  may  be  of  any  form.  The  subroutine  obtains  the  integer  array 
length  N  by  counting  the  arguments. 

Calling  Sequence — BLDAR  Y(  A( D V ),X  1  ,X2,X3, . . .  ,XN) 


Subroutine  BRKARY  or  BKARAD 

Purpose— These  subroutines  will  distribute  values  from  within  an  array  to  a  variable 
number  of  arguments  in  the  order  listed.  The  first  places  the  value  into  the  location 
while  the  second  adds  it  to  what  is  in  the  location.  Respective  operations  are 

X;  -  Aj,  i  =  1,  N 
or  X;  =  Xj  +  Aj,  i  =  1,  N. 

Restrictions— Floating-point  numbers  must  be  used  for  BKARAD.  The  integer  array 
length  N  is  obtained  by  the  routines  by  counting  the  number  of  arguments. 

Calling  Sequence— BRKARY(A(DV),X1,X2,X3 . XN)  or 

BKARAD(  A(DV),X1,X2,X3 . XN) 


Subroutine  STFSEP  or  SCALE 

Purpose— Subroutine  STFSEP  will  place  a  constant  value  into  a  variable  number  of 
locations.  Subroutine  SCALE  will  utilize  a  constant  value  to  multiply  a  variable  number 
of  arguments,  each  having  a  location  for  the  product.  The  respective  operations  are 

Xj  =  Y,  i  =  1,  2,  3 . N 

or  Xj  =  Y  *  Zj,  i  =  1,  2,  3 . N. 
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Restrictions— STFSEP  may  be  used  to  move  any  desired  value,  buc  SCALE  can  only 
be  used  for  floating-point  numbers. 

Calling  Sequence- STFSEP(Y,X1,X2,X3, . . .  ,XN)  or 

SCALE(Y,X1,Z1,X‘2,Z2 . XN,ZN) 


Subroutine  STFSEQ  or  STFSQS 

Purpose— Both  subroutines  will  stuff  a  constant  data  value  into  a  specified  length 
array  or  group  of  sequential  locations.  STFSEQ  expects  the  constant  data  value  to  be  in 
the  first  array  location,  while  STFSQS  requires  it  to  be  supplied  as  an  additional  argu¬ 
ment.  The  respective  operations  performed  are 

Ai  =  Aj,  i  =  2,  N 

or  Aj  =  B,  i  =  1,  N. 

Restrictions— N  must  be  an  integer,  but  the  constant  data  value  may  be  integer,  either 
floating  point  or  alphanumeric. 

Calling  Sequence — STFSEQ(  A( D V ),N )  or  STFSQS(B,N,A(DV)) 


Subroutine  SUMARY 

Purpose— SUMARY  is  used  to  sum  an  array  of  floating-point  values: 

S  =  2  Aj,  i  =  1,  N. 

Restrictions— The  values  to  be  summed  must  be  floating-point  numbers  and  the  array 
length  N  must  be  an  integer. 

Calling  Sequence— SUM ARY(N,A(DV),S) 


Subroutine  MAXDAR  or  MXDRAL 

Purpose— These  subroutines  will  obtain  the  absolute  maximum  difference  between 
corresponding  f  lements  of  two  arrays  of  equal  length  N.  The  array  values  must  be 
floating-point  numbers.  The  operation  performed  is 

D  =  ‘Ai  “  Bi'maxi  i=i.N- 

Subroutine  MXDRAL  also  locates  the  position  P  between  1  and  N  where  the  maximum 
occurs. 

Restrictions— The  N  argument  must  be  an  integer.  The  D  and  P  arguments  are  re¬ 
turned  as  floating-point  numbers. 

Calling  Sequence— MAXDAR(N,A(DV),B(DV).D)  or  MXDRAL(N,A(DV),B(DV),D,P) 
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Subroutine  ARYINV  or  ARINDV 

Purpose— Subroutine  ARYINV  will  invert  each  element  of  an  array  in  its  own  loca¬ 
tion.  Subroutine  ARINDV  will  divide  each  element  of  an  array  into  a  constant  value  to 
form  a  new  arra>  Their  respective  operations  are 

A;  =  1.0/ A;,  i  =  1,  N 
or  A;  =  B/C;,  i  =  1,  N. 

Restrict'ons— All  data  values  must  be  floating-point  numbers.  The  array  length  N 
must  be  an  integer. 

Calling  Sequence— ARYINV(N,A{DV))  or  ARINDV(N,C(DV),B,A(DV)) 

The  ARINDV  answer  array  may  be  overlayed  into  the  input  array  area. 


Subroutine  ADDINV  or  ADARIN 

Purpose— Subroutine  ADDINV  will  calculate  one  over  the  sum  of  the  inverses  of  a 
variable  number  of  arguments.  Subroutine  ADARIN  will  calculate  one  ovrr  the  sum  of 
inverses  of  an  array  of  values.  These  subroutines  are  useful  for  calculat  .ig  the  effective 
conductance  of  series  conductors.  Their  respective  operations  are 

Y  =  1.0/(l./X1+l./X2+...+l./XN),  N>2 
or  Y  =  l.O/Xd./Xj),  i  =  1,  N. 

Restrictions— All  data  values  must  be  floating-point  numbers.  The  array  length  N 
must  be  an  integer. 

Calling  Sequence- ADDINV<X1,X2,X3, . . .  ,XN,Y)  or  ADARIN(N,X(DV),Y) 


Subroutine  STOARY  or  ARYSTO 

Purpose— These  subroutines  will  place  a  value  into  or  take  a  value  out  of  a  specific 
array  location,  respectively.  Their  respective  operations  are 

Aj  =  X,  i  =  N,  N  >  0 

or  X  =  Aj,  i  =  N,  N  >  0. 

Restrictions— The  values  may  be  anything,  but  N  must  be  an  integer. 

Calling  Sequence- STOARY(N,X,A(DV))  or  ARYSTO(N,X,A(DV)) 
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Subroutine  SCLDEP  or  SCLIND 

Purpose— These  subroutines  will  multiply  the  dependent  or  independent  variables  of 
a  doublet  interpolation  array,  respectively.  Their  respective  operations  are 

A;  =  X  *  Aj,  i  =  3,  5,  7 . N+l 

or  Aj  =  X  *  A;,  i  =  2,  4,  6, . . . ,  N. 

Restrictions— All  values  must  be  floating  point.  The  arrays  must  contain  the  length 
integer  count  as  the  first  value,  which  must  be  even. 

Calling  Sequence- SCLDEP(A<1C),X)  or  SCLIND(A(IC),X) 


Subroutine  SLDARY  or  SLDARD 


Purpose— These  subroutines  are  useful  for  updating  fixed-length  interpolation  arrays 
during  a  transient  analysis.  The  array  data  values  are  moved  back  one  or  two  positions, 
the  first  one  or  two  values  are  discarded,  and  the  last  one  or  two  values  updated,  respec¬ 
tively.  The  “sliding  array”  thus  maintained  can  then  be  used  with  standard  interpolation 
subroutines  to  simulate  transport  delay  phenomena.  Their  respective  operations  are 


A 

and  A 
or  A 
and  A 


*  Ai+ 1  * 

-  X, 

-  Ai+2» 

=  X  and  Ai+1  =  Y, 


=  2,  N 
=  N  +  1 
=  2,  N-l 
=  N. 


Restrictions— The  addressed  arrays  must  have  the  array  integer  count  N  as  the  first 
value.  For  SLDARD,  N  must  be  even. 


Calling  Sequence— SLDARY {X,A(IC))  or  SLDARD(X,Y,A(IC)) 


Subroutine  SPLIT  or  JOIN 

Purpose— These  subroutines  separate  a  doublet  array  into  two  singlet  arrays  or  com¬ 
bine  two  singlet  arrays  into  a  doublet  array  respectively.  Their  respective  operations  are 


A2i-1  > 

i  = 

1,  N 

A2i> 

i  = 

1,  N 

=  Bi> 

i  = 

1,  N 

c«. 

i  = 

3,  N. 

Restrictions— The  arrays  may  contain  any  values,  but  N  must  be  an  integer.  N  is  the 
length  of  the  B  and  C  arrays,  and  the  A  array  must  be  of  length  2N. 

Calling  Sequence— SPLIT(N,A(DV),B(DV),C(DV))  or 
JO!N(N,B(DV),C{DV),A(DV)) 
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Purpose— This  subroutine  applies  interpolation  subroutine  D1D1DA  to  two  singlet 
arrays  to  obtain  an  array  of  dependent  variables  vs  an  array  of  independent  variables.  It 
is  extremely  useful  for  obtaining  singlet  arrays  of  various  dependent  variables  with  a  cor¬ 
responding  relationship  to  one  singlet  independent  variable  array.  The  dependent  variable 
arrays  thus  constructed  can  then  be  operated  on  by  array  manipulation  subroutines  in 
order  to  form  composite  or  complex  functions.  Doublet  arrays  can  first  be  separated 
with  subroutine  SPLIT  and  later  reformed  with  subroutine  JOIN. 


Restrictions— All  data  values  must  be  floating  point  except  N,  which  must  be  the 
integer  length  of  the  array  to  be  constructed.  The  arrays  fed  into  D1D1DA  for  interpola¬ 
tion  must  start  with  the  integer  count.  X  is  for  independent  and  Y  is  for  dependent.  I 
is  for  input  and  O  for  output. 

Calling  Sequence— SPREAD(N,X(IC),Y(IC),XI(DV),YO(DV)) 


Subroutine  QMETER  or  RDTNQS  or  QMTRI  or  QFORCE 

Purpose— These  subroutines  are  generally  used  for  calculating  flow  rates.  Their 
respective  operations  are 

A  =  B  *  (C-D) 

or  A  =  B  *  ((0+460.)“*  -  (D+460.)4) 
or  A;  =  Bj  *  (Cj-C|+  j ),  i  =  1,  N 

or  A;  =  Bj  *  (CrDj),  i  =  1,  N. 

Restrict  ions— All  values  must  be  floating  rN.N  numbers  except  the  array  length  N, 
which  must  be  an  integer. 

Calling  Sequence- QMETER(C,D,B,A)  or  RDTNQS(D,C,B,A)  or 
QMTRI(N,C(DV),B(DV),A(DV))  or 
QFORCE(N,C(DV),D(DV),B(DV),A(DV)) 


Subroutine  QINTEG  or  QINTGI 

Purpose— These  subroutines  perform  a  simple  integration.  They  are  useful  for  ob¬ 
taining  the  integrals  of  flow  rates  calculated  by  QMETER,  RDTNQS,  QMTRI,  or  QFORCE. 
Their  respective  operations  are 

S  =  S  +  Q  *  DT 
or  Sj  =  S;  +  Qj  *  DT,  i  =  1,  N. 

Restrictions— All  values  must  be  floating-point  numbers  except  N  which  must  be  an 
integer.  Control  constant  DTIMEU  should  be  used  for  the  step  size  when  doing  an  inte¬ 
gration  with  respect  to  time.  These  subroutines  should  be  called  in  Variables  2. 

Calling  Sequence— QINTEG(Q,DT,S)  or  QINTGI(N,Q(DV),DT,S(DV)) 
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Subroutine  CINSIN  or  SINARY 

Purpose  —These  subroutines  obtain  the  sine  function  of  an  angle  or  an  array  of 
angles.  Their  respective  operations  are 

A  =  sin  (B) 

or  Aj  =  sin  (B;),  i  =  1,  N. 

Restrictions— All  angles  must  be  in  radians.  AH  values  must  be  floating-point  num¬ 
bers  except  N,  which  must  be  an  integer. 

Calling  Sequence— CINSIN(B,  A)  or  SINARY(N,B(DV),A(DV)) 


Subroutine  CINCOS  or  COSARY 

Purpose— These  subroutines  obtain  the  cosine  function  of  an  angle  or  array  of 
angles.  Their  respective  operations  are 

A  =  (B) 

or  A;  =  cos  (Bj),  i  •=  1,  N. 

Restrictions— All  angles  must  be  in  radians.  All  values  must  be  floating-point  num¬ 
bers  except  the  array  length  N,  which  must  be  an  integer. 

Calling  Sequence— ClNCOS(B,  A)  or  COSARY(N,B(DV),A(DV)) 


Subroutine  CINTAN  or  TANARY 

Purpose— These  subroutines  obtain  the  tangent  function  of  an  angle  or  array  of 
angles.  Their  respective  operations  are 

A  =  tan  (B) 

or  Aj  =  tan  (B;),  i  -  1,  N. 

Restrictions— All  angles  must  be  in  radians.  Ml  values  must  be  floating  point  num¬ 
bers  except  the  array  length  N,  which  must  be  an  integer. 

Calling  Sequence- CINTAN(B.A)  or  TAN  ARY(N,B(DV),A(DV)) 


Subroutine  ARCSIN  or  ASNARY 

Purpose— These  subroutines  obtain  the  angle  corresponding  to  a  sine  function  value 
or  array  of  sine  values.  Their  respective  operations  are 

A  =  sin_,(B) 

or  Aj  =  sin-1  (Bj),  i  -  1,  N. 

Restrictions— The  angles  are  returned  in  radians  with  the  limits  A^jt/2.  All 

values  must  be  floating  point  except  for  the  array  length  N,  which  must  be  an  integer. 
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Calling  Sequence- ARCSIN(B,A)  or  ASNARY(N,B(DV),A(DV)) 

Subroutine  ARCCOS  or  ACSARY 

Purpose— These  subroutines  obtain  the  angle  correspondin':  to  a  cosine  function 
value  or  array  of  cosine  values.  Their  respective  operations  are 

A  =  cos-,(B) 

or  A,  -  cos_1(B[),  i  =  1,  N. 

Restrictions— The  angles  are  returned  in  radians  with  the  limits  0  <  A  <  n.  All 
vaiues  must  be  floating-point  numbers  except  for  the  array  length  N,  which  must  be  an 
integer. 

Calling  Sequence — ARCCOSf B, A )  or  ACSARY(N ,B(DV),A(DV») 

Subroutine  ARCTAN  or  ATNARY 

Purpose— These  subroutine  ootain  the  angle  corresponding  to  a  tangent  function 
value  of  array  of  tangent  values.  Their  respective  operations  are 

A  =  tan“!(B) 

or  A-,  =  tan-1tBj),  i  =  1.  N. 

Restrictions— The  angles  are  returned  in  radians  with  the  limits  — 7r/2  <  A  <  tt/2.  All 
values  must  be  floating-point  numbers  except  the  array  length  N,  which  must  be  an  integer. 

Calling  Sequence-  AROTAN(B.A)  or  ATNARY (N,B( D V ),A( D V ) ) 

Subroutine  EXPNTL  or  ARYEXP  or  EXPARY 

Purpose— These  subroutines  perform  an  exponential  operation.  Their  respective 
operations  are 

A  =  Bc 

or  A,  =  B,c,  I  =  1.  N 
or  Aj  =  Bf>,  1  =  1,  N. 

Restrictions— All  '.dues  must  be  positive  floating-point  numbers  except  N.  which 
must  be  an  integer 

Calling  Sequence -EXVm'UC,B. A )  or  ARYEXP(N.C.B(DV),A(DV))  or 
EXP  A  R  Y(  N  ,C(  D  V  ),B(  D  V ).  A(  D  V ) ) 
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Subroutine  LOOT  or  LOGTAR 

Purpose  -These  subroutines  obtain  the  base  10  log  function  of  a  number  or  array  of 
numbers.  Then  respective  opt  rations  are 

A  =  logj  0(  B) 

or  A,  =  log, 0( B, ).  i  =  1.  N. 

Restrictions  —  All  values  must  be  positive  floating-point  numbers  except  N,  which 
must  be  an  integer. 

Calling  Sequence  -LOGTtB.Ai  or  LOGTAR(N.B(r  V).A(DV)) 

Subroutine  LOCK  or  LOGEAR 

Purpose  -  These  subroutines  obtain  the  base  e  log  function  of  a  number  or  array  of 
numbers  Their  respective  operations  are 

A  =  lcge(Bi 

or  A,  =  log(.(B, ).  l  -  1.  N. 

Restrictions  -  All  values  must  be  positive  floating-point  numbers  except  N,  which 
must  be  an  int<*ger. 

Calling  Sequence  —  LOGK( B.A  I  or  LOGEAR(N.B( DV ),A( DV ) ) 

Subroutine  SQROOT  or  SQROTI 

Purpose  These  subroutines  obtain  the  square  root  of  a  number  or  array  of  numbers, 
respectively.  Their  respective  operations  are 

A  =  ->-s/  o 

or  A,  ~  +  \/B,.  i  1.  N. 

Restrictions— ‘the  A  and  B  values  must  be  floating-point  numbers  The  N  must  be 
an  integer. 

Calling  Sequence  SQROOTt B.Ai  or  SQKOTIlN.B<r>V|.A<!)Vn 
Subroutine  CMPXSR  or  CSQRI 

Purpose  -These  subroutines  obtain  the  complex  square  root  of  a  complex  number 
or  an  array  of  complex  numbers,  respectively.  Their  respective  operations  are 

A  +  iB  -  y  (  4  >D.  i  "  \'  - 1 

or  A,  ♦  iB,  =  \  j  *  >!>,.  j  --  1.  V 
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Restrictions— All  numbers  must  be  floating-point  except  N,  which  must  be  an 
integer. 

Calling  Sequence— CMPXSR(C,D, A, B)  or  CSQRI(N,C(DV),D(DV),A(DV),B(DV)) 

Subroutine  CMPXMP  or  CMPYI 

Purpose— These  subroutines  will  multiply  two  complex  numbers  or  the  correspond¬ 
ing  elements  of  arrays  of  complex  numbers.  Their  respective  operations  are 

A  +  iB  =  (C  +  iD)*(E  +  iF),  i  =  v'-T 

or  Aj  +  iB,  =  (Cj  +  iDj)^(Ej  +  iFj),  j  =  1,  N 

Restrictions— All  numbers  must  be  floating  point  except  for  N,  which  must  be  an 
integer. 

Calling  Sequence— CMPXMP(C,D,E,F,A,B)  or 

CMPY1(N,C(DV).D(DV),E(DV),F(DV),A(DV),B(DV)) 

Subroutine  CMPXDV  or  CDIVI 

Purpose— These  subroutines  will  divide  two  complex  numbers  or  the  corresponding 
elements  of  arrays  of  complex  numbers.  Their  respective  operations  are 

A  +  iB  =  (C  +  iD)/(E  +  iF),  i  =  ^ 

or  Aj  +  iBj  =  (Cj  +  iDj)/(E,  +  iF,).  j  =  l.N 

Restrictions— All  numbers  must  be  floating  point  except  for  N,  which  must  be  an 
integer. 

Calling  Sequence- CMPXDV(C,D,E,F,A,B)  or 

CDIVI(N’,C(DV>,D(DV  ),E(DV),F(DV),A(DV),B(DV)) 

Subroutine  NEWTRT  or  N^WRT  l 

Purpose— These  subroutines  utilize  Newton’s  method  to  obtain  one  root  of  a  cubic 
or  quartic  equation,  respectively.  The  root  must  be  in  the  neighborhood  of  the  supplied 
initial  guess,  and  up  to  100  iterations  are  performed  in  order  to  obtain  an  answer  witb'n 
the  specified  tolerance.  If  the  tolerance  is  not  met,  an  answer  of  1G38  is  returned.  The 
respective  equations  are 

f(X)  =  A 1  +  A  2  *  X + A  3  *  X  2 + A4  *  X  3  =  O.OiT 
or  g( X )  =  Al+A2+X  * .\3*X2+A-1*X3+A5*X'1  =  O.OiT 

wI/tc  X  starts  as  the  initial  guess  R!  and  finishes  as  the  final  answer  RF.  T  is  the 
toleience. 
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Restrictions— All  data  values  must  he  floating-point  numbers. 

Calling  Sequence— NE\VTRT(A(DV),T,RI,RF)  or  NEWRT4(A(DV),T,RI,RF) 


Subroutine  PLYNML  or  PLY  ARY 

Purpose— These  subroutines  calculate  Y  from  the  following  polynomial  equation: 

Y  =  A1+A2*X+A3*X2+A4*X3+  ...  +AN*XN_1. 

The  number  of  terms  is  variable,  but  all  the  A  coefficients  must  be  entered  no  matter 
what  their  value. 

Restrictions— All  values  must  be  floating-point  numbers  except  the  number  of  coeffi¬ 
cients  N,  which  must  be  an  integer. 

Calling  Sequence — PL YNM L( X, A 1 , A2, A3 . AN,Y)  or  PLYARY(N,X,A(DV),Y) 


Subroutine  SMPINT  or  TRPZD 

Purpose— These  subroutines  perform  area  integrations  by  Simpson’s  rule  and  the 
trapezoidal  rule,  respectively.  Simpson’s  rule  requires  that  an  odd  number  of  points  be 
supplied.  If  an  even  number  of  points  is  supplied,  SMPINT  will  apply  the  trapezoidal  rule 
to  the  last  incremental  area  but  Simpson’s  rule  elsewhere.  The  respective  operations  are 

A  =  DX*( Yl+4 Y2+2Y3+4Y4+  . . .  +YN)/3 
or  A  =  DX*(Ylf2Y2+2Y3+2Y4+...+YN)/2. 

Restrictions— The  DX  increment  must  be  uniform  between  all  the  Y  points.  All 
values  must  be  floating  point  except  N,  which  must  be  an  integer. 

Calling  Seq uence — SM PI NT( N , DX , Y ( D V ), A )  or  TRPZD(N,DX,Y(DV),A) 


Subroutine  TRPZDA 

Purpose  —  This  subroutine  performs  area  integietion  by  the  trapezoidal  rule.  It 
should  be  used  where  the  DX  increment  is  not  uniform  between  the  Y  values  but  tne 
corresponding  X  value  for  each  Y  value  is  known.  The  operation  performed  is  ,»s  follows: 

A  =  ji:(Xi-X,-l)*(VYl-l),  i  =  2.  N. 

Restrictions— All  values  must  be  floating-point  numbers  except  the  army  length  N, 
which  must  be  an  integer. 

Calling  Sequence  -TRPZDA(N.X(DV).Y(DV),A) 
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Subroutine  PRESS  or  SPRESS 

Purpose— These  routines  are  useful  for  impressing  nodal  pressures  in  one-dimensional 
flow  paths  once  the  entry  pressure  PI,  path  conductance  G,  and  flow  rate  W  are  known. 
The  respective  equations  are 

P2  =  Pl-W/G 

or  Pli+1  =  Plj-W/G,,  i  =  1,  2,  3 . N. 

Restrictions— F or  SPRESS,  the  pressures  and  conductors  must  be  sequential  and  in 
ascending  order;  the  number  of  pressure  points  to  be  calculated  must  be  supplied  as  the 
integer  N. 

Calling  Sequence — PR ESS( PI -  W ,G ,P2 )  or  SPRESS(N,P1(DV),W,G(DV)) 


Subroutine  EFFG 

Purpose— Subroutine  EFFG  is  a  pressure  network  of  the  type  in  Fig.  11. 


Where  the  values  of  tne  identified  elements  are  known,  this  subroutine  will  calculate  the 
effective  conductance  GE  from  PI  to  P2.  Any  interconnections  may  occur  in  the  space, 
but  only  P2,  P3  and  P4  may  be  on  the  boundary  and  no  elements  may  cross  it.  The 
equation  utilized  is 

GE  =  (G1  *(P1-P3)  +  G2*(P]  -P4))/(P1-P2). 

Restrictions-  See  above.  May  not  be  used  where  capacitors  appear  on  the  internal 
nodes. 

Calling  Sequence  -EFFG(P1.P2,P3,P4.G1.G2,GE) 


Subroutine  EFFE.MS 

Purpose  -Tin.  subroutine  calculates  the  effective  emissivity  E  betw^n  parallel  flat 
plates  by  the  following  equation: 

V  =  1. 0/(1. 0/El  +  1.0/E2  -  1.0), 

where  El  and  E2  are  the  ennssivities  of  the  two  surfaces  under  consideration. 
Restrictions  -  Argumi  .Us  must  be  floatmg-poin*  numbers. 


Calling  Sequence  EF FEMS(  K I  ,E2.K ; 
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Output  Subroutines 


Name  Page 

STNDRD,  PRNTMP,  PRINT,  PRINTL  .  83 

PRINTA,  PRNTMA,  PUNCHA  .  84 

TPRINT,  READ,  WRITE,  EOF,  REWIND .  85 


Subrou  ie  STNDRD  or  PRNTMP 

Purpose— Subroutine  STNDRD  causes  a  line  of  output  to  be  printed  giving  the 
present  time,  the  last  time  step  used,  the  most  recent  CSGMIN  value,  the  maximum  dif¬ 
fusion  temperature  change  calculated  over  the  last  time  step,  and  the  maximum  relaxation 
change  calculated  over  the  last  iteration.  RNN  refers  to  the  relative  node  number  on 
which  something  occurred.  The  line  of  output  looks  as  follows: 

***** 

TIME _ DTIMEU _ CSGM  JN(  RNN ) _ DTMPCC’t  RNN ) _ ARLXCC(RNN ) _ 

Subroutine  PRNTMP  internally  calls  on  STNDRD  and  also  lists  the  temperature  of  every 
node  in  the  network  according  to  relative  node  number.  The  relative  node  number  vs 
actual  n<  de  number  dictionary  printed  out  with  the  input  data  should  he  consulted  to 
determine  temperature  locations  on  the  thermal  network  model. 

Restrictions— No  arguments  are  required  or  allowed.  These  subroutines  should  be 
used  with  network  problems  only. 

Calling  Sequence- STNDRD  or  PRNTMP 


Subroutine  PRINT  or  PRINTL 

Purpose  —  These  subroutines  allow  individual  floating-point  numbers  to  be  printed  out. 
The  arguments  may  reference  tem|>eraturo,  capacitance,  source  locations,  conductors, 
constants,  or  unique  array  locations.  In  addition,  subroutine  PRINTL  allows  each  value 
to  he  preceded  or  labeled  by  a  6-character  alphanumeric  word.  The  number  of  arguments 
is  variable,  but  the  “label”  array  used  for  PRINTL  should  contain  a  label  for  each 
argument. 
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Restrictions — These  subroutines  do  not  call  on  STNDRD.  The  user  may  call  on  it  if 
he  desires  time  control  information.  Any  control  constant  may  be  addressed  in  order  to 
see  what  its  value  is;  integers  must  first  be  floated. 

Calling  Sequence  -PR  1NT(T,(\Q.G.K . A+)  or  PRINTL(LA(DV),T,C,Q,G,K,...,A+) 


Subroutine  PRINTA 

Purpose— 'Phis  subroutine  allows  the  user  to  print  out-  an  array  of  values,  five  to  the 
line.  The  integer  array  length  N  and  tin*  first  data  value  location  must  be  specified.  Each 
value  receives  an  indexed  label.  The  user  must  supply  a  6-cb  jracter  alphanumeric  word  L 
to  be  used  as  a  common  label  and  an  integer  value  M  to  begin  the  index  count. 

Restrictions— The  array  values  to  be  printed  must  be  floating-point  numbers. 

Calling  Sequence  -  PR  1  N'T  A  ( L.  A(  DV  ),N,M ) 

If  the  label  was  the  work  TEMP,  N  was  3,  and  M  was  6,  the  line  of  output  will  look  as 
follows: 

TEMP(  (>)  valueTKMP  (  7 lvalue  TEMP  (  8 lvalue 


Subroutine  PRN’TMA 

Purpose— This  subroutine  allows  the  user  to  print  out  up  to  10  arrays  in  a  column 
format.  The  individual  elements  are  not  labeled,  but  each  column  receives  a  2-line  head¬ 
ing  of  12  alphanumeric  characters  each.  The  2-line  heading  must  be  supplied  as  a  single 
array  of  four  words,  six  characters  each.  The  user  must  supply  the  starting  location  of 
each  label  array  and  value  array.  The  number  of  values  in  each  value  array  must  agree 
and  be  supplied  as  the  integer  N.  The  value  arrays  must  contain  floating-point  numbers. 

Restrictions  —  Labels  must  be  alphanumeric,  while  values  must  be  floating  point.  All 
floating-point-value  arrays  must  contain  the  same  number  of  values. 

Calling  .Ser/f/ence— PRNTMA(N,LA1  (DV ),VA1  (DV ),LA2(DV ),VA2(DV), . . .) 


Subroutine  PUNCH  A 

Purpose -This  subroutine  enables  a  user  to  punch  out  an  array  of  data  values  in  any 
desired  format.  The  F  argument  must  reference  a  FORTRAN  format  which  has  been 
input  as  an  array,  including  the  outer  parentheses  but  deleting  the  word  format.  The 
second  argument  must  address  the  first  data  value  of  the  array  of  sequential  values.  The 
third  argument  N  must  be  the  integer  number  of  data  values  in  the  array. 

Restrictions  -Punched  cards  must  lie  asked  for  on  the  job  request  form. 

Calling  Sequence  —  PUNCH A(  F(  DV  ),A(  DV  ),N ) 
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Subroutine  TPRINT 

Purpose— Subroutine  TPRINT  makes  a  call  to  STNDRD,  then  lists  the  actual  node 
number  and  corresponding  temperature  for  every  node  in  a  network. 

Restrictions— 'This  subroutine  may  be  called  from  any  of  the  operations  blocks. 

Calling  Sequence— TPRINT 


Subroutine  READ  or  WRITE 

Purpose— These  subroutines  enable  the  user  to  read  and  write  array,  of  data  as  binary 
information  on  magnetic  tape.  The  first  argument  L  must  be  the  integer  number  of  the 
logical  tape  being  addressed,  The  second  argument  X  must  address  the  first  data  value 
ol  the  array  to  be  written  out  or  the  starting  location  for  data  to  be  read  into.  The  third 
argument  N  must  be  an  integer.  For  WRITE  it  is  the  number  of  data  values  to  be  written 
on  tape  as  a  record.  For  READ  it  is  the  number  of  data  varies  to  be  read  in  from  tape 
from  the  next  record,  not  necessarily  the  entire  record. 

Restrictio ns— The  user  should  check  section  VII  to  determine  which  logical  units 
are  available  and  control  card  requirements.  All  processed  information  must  be  in  binary. 

Calling  Sequence— READ(L,X(DV),N)  or  WRITE(L,X(DV),N) 


Subroutine  EOF  or  REWIND 

Purpose— These  subroutines  enable  the  user  to  write  end  of  file  marks  on  magnetic 
tape  and  to  rewind  them.  They  are  generally  used  in  conjunction  with  subroutines  READ 
and  WRITE  discussed  above.  The  single  argument  L  must  be  the  integer  logical  tape 
number  of  the  unit  being  cctiviated. 

Restrictions— The  user  should  check  section  VII  to  determine  available  logical  units. 

Calling  Sequence -EOF  (L)  or  REWIND  (L) 
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Matrix  Subroutines 


Name  Page 

ZERO,  ONES,  UNITY,  SIGMA,  GENALP,  GENCOL .  87 

SHIFT,  REFLCT,  SHUFL,  COLMAX,  COLMIN  .  88 

ELEADD,  ELESUB,  ELEMUL,  ELEDIV,  ELEINV, 

i'.FSIN,  EFASN  .  89 

EFCOS,  EFACS,  EFTAN,  EFATN,  EFLOG,  EFSQR .  90 

EFEXP,  EFPOW,  MATRIX,  SCALAR,  DISAS,  ASSMBL .  91 

DIAG,  COLMLT,  ROWMLT,  ADDALP,  ALPHAA,  AABB  .  92 

BTAB,  INVRSE,  MULT .  93 

TRANS,  POLMLT,  POLVAL,  PLYEVL .  94 

POLSOV,  JACOBI,  MODES  .  95 

MASS  .  96 

STIFF,  LIST  .  97 

PUNCH,  Matrix  Data  Storage  and  Retrieval,  CALL,  FILE, 

ENDMOP,  LSTAPE  .  98 


NOTE:  All  of  the  ab« -ubroutines  require  that  matrixes  be  entered  as  positive  num¬ 
bered  arrays  having  the  ..ceger  number  of  rows  and  columns  as  the  first  two  data  values 
followed  by  the  floating-point  element  values  in  row  order.  The  above  package  of  sub¬ 
routines  is  referred  to  as  MOPAS,  for  Matrix  Oriented  Production  Assembly  System. 
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Subroutine  ZERO  or  ONES 

Purpose— These  subroutines  generate  a  matrix  [Z]  such  that  every  element  is  zero  or 
one,  respectively. 

Restrictions— 'The  matrix  to  be  generated  mmt  contain  exactly  enough  space  in  addi¬ 
tion  to  having  the  integer  number  of  rows  and  columns  as  the  first  two  data  values.  The 
NR  and  NC  arguments  are  the  integer  number  of  rows  and  columns,  respectively. 

Calling  Sequence— ZERO(NR,NC,Z(IC))  or  ONES(NR,NC,Z(IC)) 


Subroutine  UNITY  or  SIGMA 

Purpose— These  are  square  matrix  generation  subroutines.  UNITY  generates  a  square 
matrix  such  that  the  main  diagonal  elements  are  one  and  all  other  elements  are  zero. 
SIGMA  generates  a  square  matrix  such  that  all  elements  on  and  below  the  main  diagonal 
are  one  and  the  remaining  elements  are  zero. 

Calling  Sequence— UNITY(N,Z(IC))  or  S!GMA(N,Z(IC)) 

Restrictions— The  matrix  [Zj  to  be  generated  must  contain  exactly  enough  space  in 
addition  to  having  the  integer  number  of  rows  and  columns  as  the  first  two  data  values. 
The  integer  numl>er  of  rows  and  columns  are  equal  and  must  be  input  as  the  argument  N. 


Subroutine  GEN  ALP  or  GENCOL 

Purpose— These  are  special  matrix  generation  subroutines.  GENALP  will  generate  a 
matrix  such  that  every  element  is  equal  to  a  constant  C.  GENCOL  will  generate  a  column 
matrix  such  that  the  first  element  is  equal  to  XI  and  the  last  element  is  equal  to  X2.  The 
intermediate  elements  receive  equally  incremented  values  such  that  a  linear  relationship  is 
established  between  row  number  and  element  value. 


Restrictions— The  NR  and  NC  arguments  refer  to  the  integer  number  of  rows  and 
columns,  respectively.  XI,  X2,  and  C  must  be  floating-point  values.  The  generated 
matrixes  must  contain  exactly  enough  space  in  addition  to  having  the  integer  number  of 
rows  and  columns  as  the  first  two  data  values. 

Calling  Sequence— GENALP(NR,NC.C,Z( IC ) )  or  SENCOL( X 1  ,X2,N R,7.(  1C ) ) 
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Subroutine  SHIFT  or  REFLCT 

Purpose— These  subroutines  may  be  used  to  move  an  entire  matrix  from  one  location 
to  another.  SHIFT  moves  the  matrix  exactly  as  is  and  REFLCT  moves  it  and  reverses 
the  order  of  the  elements  within  each  column.  The  last  element  in  each  column  becomes 
the  first  and  the  first  becomes  the  last,  etc. 

REFLCT  uses  three  dynamic  storage  locations  plus  an  additional  one  for  each  row. 

Restrictions— The  matrixes  must  be  of  identical  size,  and  the  integer  number  of  rows 
and  columns  must  be  the  first  two  data  values.  The  |Z]  matrix  may  be  overlayed  into 
the  [AJ  matrix. 

Calling  Sequence —SHI  FT( A(  1C j,Z( IC) )  or  REFLCT(A(1C),Z(IC)) 


Subroutine  SHUFL 


Purpose— This  subroutine  allows  the  user  to  reorder  the  size  of  a  matrix  as  long  as 
the  total  number  of  elements  remains  unchanged.  The  row  order  input  matrix  [A]  is 
transposed  to  achieve  column  order  and  then  reformed  as  a  vector  by  sequencing  the 
columns  in  ascending  order.  This  vector  is  then  reformed  into  a  column  order  matrix  by 
taking  a  column  at  a  time  sequentially  from  the  vector.  The  newly  formed  column  ma¬ 
trix  is  then  transposed  and  output  as  the  row  order  matrix  [ Z] . 

Restrictions— The  matrixes  must  be  identical  in  size  and  have  their  respective  integer 
number  of  rows  and  columns  as  the  first  two  data  values.  The  number  of  rows  time 
columns  for  |  A]  must  equal  the  number  of  rows  times  columns  of  [  Z] . 

Calling  Sequence— SHUFL(A(IC),Z(1C)) 


Subroutine  COLMAX  or  COLMIN 


Purpose  --These  subroutines  search  an  input  matrix  to  obtain  the  maximum  or  mini¬ 
mum  values  within  each  column,  respectively.  These  values  are  output  as  a  single  row 
matrix  |Z|  having  as  many  columns  as  the  input  matrix  [Aj. 

Restrictions— Each  matrix  must  have  its  integer  number  of  rows  and  columns  as  the 
first  two  data  values. 

Calling  .Sequence— COLMAX(A(IC).Z(JC)>  or  COLMIN(A(IC),Z(JC)) 


88 


NRL  REPORT  7656 


Subroutine  ELEADD  or  ELESUB 

Purpose— These  subroutine*:  add  or  subtract  the  corresponding  elements  of  two 
matrixes,  respectively; 

m*n  m*n  m*n 

[Z)  =  [A]  ±  [B],  z;j  =  ajj  ±  bn. 

Restrictions— All  matrixes  must  be  of  identical  size  and  have  the  integer  number  of 
rows  and  columns  as  the  first  two  data  values.  The  [ZJ  matrix  may  be  overlayed  into 
the  [A]  or  |B]  matrix. 

Calling  Sequence— ELEADD(A(IC),B(1C),Z(IC))  or  ELESUB(A(1C),B(1C),Z(IC)) 


Subroutine  ELEMUL  or  ELEDIV 

Purpose— These  subroutines  multiply  or  divide  the  corresponding  elements  of  two 
matrixes,  respectively; 

m*n  m*n  m*n 
[Z]  =  [A]  */  [B],  Zjj  -  ajj  */  bjj. 

Restrictions— All  matrixes  must  be  of  identical  size  and  have  the  integer  number  of 
rows  and  columns  as  the  first  two  data  values.  The  (Z)  matrix  may  be  overlayed  into  the 
| A]  or  (B)  matrix. 

Calling  Sequence —  ELEM  U  L{  A(  IC  ),B(  IC  ),Z(  1C ) )  or  ELED1V(A(IC),B(IC),Z(JC)) 
Subroutine  ELEINV 

Purpose— This  subroutine  obtains  the  reciprocal  of  each  element  of  the  [A]  matrix 
and  places  it  in  the  corresponding  element  location  of  the  [Z]  matrix; 

Zjj  =  1.0/a,,. 

Restrictions— The  matrixes  must  be  of  identical  size  and  have  the  integer  number  of 
rows  and  columns  as  the  first  two  data  values.  The  |Z]  matrix  may  be  overlayed  into  the 
(A)  matrix. 

Calling  Sequence  —  ELE1N V(  A( IC),Z( IC)) 


Subroutine  EFSIN  or  EFASN 

Purpose— These  subroutines  perform  elementary  functions  on  all  of  the  (A)  matrix 
elements  as  follows: 

z,j  =  sin  ( atJ )  or  z(J  =  art-sin  (ajj). 
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Restrictions— The  matrixes  must  be  identical  in  size  and  have  the  integer  number  of 
rows  and  columns  as  the  first  two  data  values.  The  [Z]  matrix  may  be  overlayed  into  the 
(Aj  matrix. 

Calling  Sequence- EFSIN(A(IC),Z(IC))  or  EFASN(A(1C),Z(IC)) 


Subroutine  EFCOS  or  EFACS 

Purpose— These  subroutines  perform  elementary  functions  on  all  of  the  [A]  matrix 
elements  as  follows: 


Zjj  =  cos  (a5j )  or  zi}  =  arccos  (ajj ). 

Restrictions— The  matrixes  must  be  identical  in  size  and  have  the  integer  number  of 
rows  and  columns  as  the  first  two  data  values.  The  [Zj  matrix  may  be  overlayed  into 
the  [A]  matrix. 

Calling  Sequence- EFCOS(A(IC),Z(IC))  or  EFACS(A(IC),Z(IC)) 


Subroutine  EFTAN  or  EFATN 

Purpose— These  subroutines  perform  elementary  functions  on  all  of  the  (AJ  matrix 
elements  as  follows: 


z^  =  tan  (ajj)  or  Zjj  =  arctan  (ajji. 

Restrictions  —The  matrixes  must  be  of  identical  size  and  have  the  integer  number  of 
rows  and  columns  as  the  first  two  data  values.  The  [ZJ  matrix  may  be  overlayed  into 
the  [ A 1  matrix. 

Calling  Sequence — EFTAN'( A( IC ),Z( IC) )  or  EFATN(A(IC),ZMC)) 


Subroutine  EFLOG  or  EFSQR 

Purpose— These  subroutines  jierform  elementary  functions  on  all  of  the  [A)  matrix 
elements  as  follows: 


zu =  1op0<aij>  or  zij=V*iT 

Res  trie  tit  is— The  matrixes  must  lie  identical  in  size  and  have  the  integer  number  of 
rows  anti  columns  as  the  first  two  data  values.  All  elements  in  the  [A]  matrix  must  be 
positive. 

Calling  Sequence  —  EF LOG( A( IC ),?■( IC ) )  or  EFSQR(A(1C),Z(IC)) 
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Subroutine  EFEXP  or  EFPOW 


Purpose— These  subroutines  perform  elementary  functions  on  all  of  the  [A]  matrix 
elements  as  follows: 

Zij  =  eaij  or  zij  =  a-  ■ 

Restrictions— The  matrixes  must  be  identical  in  size  and  have  the  integer  number  of 
rows  and  columns  as  the  first  two  data  values.  The  [Z]  matrix  may  be  overlayed  into 
the  [  A 1  matrix.  The  exponent  a  may  be  an  integer  or  fion  ing-point  number.  However, 
if  any  elements  in  (A]  are  negative  then  a  must  be  an  integer. 

Calling  Sequence -EFEXP(A(IC),Z(IC))  or  EFPOW(A(IC),o,Z(IC)) 


Subroutine  MATRIX  or  SCALAR 

Purpose— Subroutine  MATRIX  allows  constant  to  replace  a  specific  matrix  element, 
and  subroutine  SCALAR  allows  a  specific  matrix  element  to  be  placed  into  a  constant 
location.  The  integers  1  and  J  designate  the  row  and  column  position  of  the  specific 
element; 


Zj;  =  C  or  C  =  Zjj. 

Restrictions— The  matrix  must  have  the  integer  number  of  rows  and  columns  as  the 
first  two  data  values.  Checks  are  made  to  insure  that  the  identified  element  is  v/ithin  the 
matrix  boundaries. 

Calling  Sequence MIC))  or  SCALAR(Z(IC),I,J,C) 


Subroutine  DISAS  or  ASSMBL 

Purpose —These  subroutines  allow  a  user  to  operate  on  matrixes  in  a  partitioned 
manner  by  disassembling  a  submatrix  [Z]  from  a  parent  matrix  (A]  cr  assembling  a  sub¬ 
matrix  1Z]  into  a  parent  matrix  j  A]  - 

Restrictions— The  I  and  J  arguments  are  integers  which  identify  (by  row  and  column 
number,  respectively)  the  upper  left-hand  corner  position  of  the  submatrix  within  the 
parent  matrix.  All  matrixes  must  have  exactly  enough  space  and  contain  the  integer 
number  of  rows  and  columns  as  the  first  two  data  values.  The  NR  and  NC  arguments  are 
the  integer  number  of  rows  and  columns.  res|>ective!y,  of  the  disassembled  submatrix.  If 
the  submatrix  exceeds  the  bounds  of  the  parent  matrix  an  appropriate  error  message  is 
written  and  the  program  terminated. 

Calling  Sequence — DIS AS( A( IC),I ,J,NR ,NC,Z(  1C) )  or  ASSMBL(Z(IC),I,J,A(IC)) 
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Subroutine  DIAG 

Purpose— Given  a  1*N  or  N*1  matrix  [V],  this  subroutine  forms  a  full  square  N*N 
matrix  ( Z] .  The  [V]  values  are  placed  sequentially  on  the  main  diagonal  of  [Z]  and  ali 
off-diagonal  elements  are  set  to  zero. 

Restrictions— Both  matrixes  must  have  exactly  enough  space  and  contain  their  integer 
number  of  rows  and  columns  as  the  first  two  data  values. 

Calling  Sequence— D1AG(V(IC),Z(IC)) 

Subroutine  COLMLT  or  ROWMLT 

Purpose— To  multiply  each  element  in  a  column  or  row  of  matrix  [A]  by  its  cor¬ 
responding  element  from  the  matrix  [V]  which  is  conceptually  a  diagonal  matrix  but 
stored  as  a  vector;  i.e.,  1*N  or  N*1  matrix.  The  matrix  [Z]  is  the  product. 

Restrictions— The  matrixes  must  have  exactly  enough  space  and  contain  the  integer 
number  of  rows  and  columns  as  the  first  two  data  values.  The  matrixes  being  multiplied 
must  be  conformable. 

Calling  Sequence— COLMLT(A(IC),V(IC),Z(IC))  or  ROWMLT(V(IC),A(IC),Z(IC)) 


Subroutine  ADDALP  or  ALPHAA 

Purpose— These  subroutines  add  a  constant  to  or  multiply  a  constant  times  every 
element  in  a  matrix; 

zsj  =  C+ajj  or  Zij  =  C*aij. 

Resttictia'.s— The  matrixes  must  have  exactly  enough  space  and  contain  the  integer 
number  of  rows  and  columns  as  the  first  two  data  values.  C  and  all  elements  must  be 
floating-point  numbers.  The  [Z]  matrix  may  be  overlayed  into  the  [A]  matrix. 

Calling  Sequence- ADDALP(C,A(IC),Z(IC))  or  ALPHAA(C,A(IC),Z(IC)) 


Subroutine  AABB 

Purpose— ‘To  sum  two  scaled  matrixes; 

m*"n  m*n  m*n 

IZ]  =  Cl  (A]  +  C2  [B]  and  zV]  =  Cl*'dV)  +  C2*bjj. 

Restrictions— All  matrixes  must  be  of  identical  size,  contain  exactly  enough  space, 
and  contain  the  integer  number  of  rows  and  columns  as  the  first  two  data  values.  The 
output  matrix  [Z]  may  be  overlayed  into  either  of  the  input  matrixes. 

Calling  Sequence- AABB(C1,A(IC),C2,B(1C),Z(1C)) 
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Subroutine  BTAB 

Purpose— To  perform  the  following  matrix  operation: 

n*n  n*m  m*m  m*n 
[ZJ  =  [B] 1  [A]  [B]  ' 

Restrictions— The  matrixes  must  be  conformable,  contain  exactly  enough  space,  and 
contain  the  integer  number  of  rows  and  columns  as  the  first  two  data  values.  Subrou¬ 
tines  MULT  am’  TRANS  are  called  on. 

This  subroutine  (due  to  MULT  and  TRANS)  uses  2*m*n+6  dynamic  storage 
locations. 

Calling  Sequence- BTAB(A(IC),B(1C),Z(1C)) 


Subroutine  INVRSE 

Purpose— To  invert  a  square  matrix; 

n*n  n*n  n*n 
given  [A],  (Z)  =  [A]'1. 

This  subroutine  requires  n  dynamic  storage  locations. 

Restrictions— The  matrixes  must  be  square,  identical  in  size,  and  contain  the  integer 
number  of  rows  and  columns  as  the  first  two  data  values.  The  output  matrix  [Z]  may 
be  overlayed  into  the  [A]  matrix. 

Calling  Sequence— INVRSE(A(IC),Z(IC)) 


Subroutine  MULT 

Purpose— To  multiply  two  conformable  matrixes  together; 

m*n  m*p  p*n 

(ZJ  =  l A 1  IB],  Zjj  =  aik*bkj. 

This  subroutine  requires  n*m  dynamic  storage  locations. 

Restrictions— The  matrixes  must  have  exactly  enough  space  and  contain  their  integer 
number  of  rows  and  columns  as  the  first  two  data  values.  If  (A]  and  [B]  are  square, 

[Z]  may  be  overlayed  into  either  of  them. 

Calling  Sequence— MULT(A(IC),B(K'),Z(IC)) 
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Subroutine  TRANS 

m*n  n*m 

Purpose— Given  a  matrix  [A],  form  its  transpose  as  [Z] . 

This  subroutine  requires  n*m  dynamic  storage  locations. 

Restrictions— Both  matrixes  must  have  exactly  enough  space  and  contain  their  inte¬ 
ger  number  of  rows  and  columns  as  the  first  two  data  values.  The  output  matrix  [Z] 
may  be  overlayed  into  the  [A]  matrix. 

Calling  Sequence— TRANS(A(IC),Z(IC)) 


Subroutine  POLMLT 

Purpose— This  subroutine  performs  the  multiplication  of  a  given  number  of  nth 
order  polynomial  coefficients  by  a  similar  number  o'  mth  order  polynomial  coefficients. 
The  polynomials  must  be  input  as  matrixes  with  the  number  of  rows  equal,  and  each  row 
receives  the  following  operation: 

(Cj  ,c2 ,c3 , . . .  ,ck )  =  (aa  ,a2 , . . .  ,an  )*(bj  ,b2 , . . .  ,b  k  =  m+n-1. 

Restrictions— The  matrixes  must  have  exactly  enough  space  and  contain  their  inte¬ 
ger  number  of  rows  and  columns  as  the  first  two  data  values. 

Calling  Sequence— POLMLT(A(IC),B(IC),((IC)) 


Subroutine  POLVAL 

Purpose— Given  a  set  of  polynomial  coefficients  as  the  first  row  of  matrix  [A],  this 
subroutine  evaluates  the  polynomial  for  the  input  complex  number  X+iY.  The  answer  is 
returned  as  U+iV. 

Restrictions— ( A]  may  be  m*n,  but  only  the  first  row  is  evaluated. 

Calling  Sequence— POLVAL(A(IC),X,Y,U,V) 


Subroutine  PLYEVL 

Purpose— Given  a  matrix  [A]  containing  an  arbitrary  number  NRA  of  nth  order 
polynomial  coefficients  and  a  column  matrix  [XJ  containing  an  arbitrary  number  NRX 
of  x  values,  this  subroutine  evaluates  each  polynomial  for  each  x  value.  The  answers  are 
output  as  a  matrix  [Z]  of  size  NRX*NRA.  Each  set  of  polynomial  coefficients  in  [ A 1 
is  a  row  in  ascending  order.  An  x  value  evaluated  for  the  polynomials  creates  a  row  in 
[Z]  where  the  column  number  agrees  with  the  polynomial  row  number. 
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Restrictions— The  matrixes  must  have  exactly  enough  space  and  contain  their  integer 
number  of  rows  and  columns  as  the  first  two  data  values. 

Calling  Sequence- PLYEVL(A(IC),Z(1C),Z(IC)) 


Subroutine  POLSOV 

Purpose— Given  a  set  of  polynomial  coefficients  as  the  first  row  in  matrix  [A], 
size  (m,n+l),  this  subroutine  calculates  the  complex  roots  which  are  returned  as  matrix 
[Z],  size  (n,2).  Column  1  contains  the  real  part  and  column  2  the  imaginary  part  of 
the  roots. 

Restrictions— This  subroutine  presently  is  limited  to  n  =  20.  It  internally  calls  on 
RTPOLY  and  utilizes  some  double  precision. 

Calling  Sequence — PO LSO V ( A( IC ) ,Z( IC ) ) 


Subroutine  JACOBI 

Purpose— This  subroutine  will  find  the  eigenvalues  [E]  and  eigenvector  matrix  [Z] 
associated  with  an  input  matrix  [A]; 

n*n  n*n  n*n  n*l 
[A]  [Z1  =  (Z|  [E]. 

This  subroutine  requires  2*n*n+6  dynamic  storage  locations. 

Restrictions— The  matrixes  must  have  exactly  enough  space  and  contain  their  integer 
number  of  rows  and  columns  as  the  first  two  data  values. 

Calling  Sequence— JACOBI(A(IC),E(IC),Z(!C)) 


Subroutine  MODES 

Purpose— This  subroutine  solves  the  dynamic  vibration  equation 

n*n  n*n  n*n  n*n  n*l 
(A)  (Z]  =  [B]  [ Z 1  m, 

W2 

where  (A]  is  the  input  inertia  matrix  associated  with  the  kinetic  energy  and  [B)  is  the 
input  stiffness  matnx  associated  with  the  strain  energy.  |Z]  is  the  output  eigenvector 
matrix  associated  with  the  frequencies  of  vibration  Wj  which  are  output  in  rad/sec  as 
[R]  and  in  hertzes  as  ( C] ;  both  [R]  and  (C]  are  n*l  matrixes. 

This  subroutine  requires  3*n*n+9  dynamic  storage  locations. 
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Restrictions— The  matrixes  must  have  exactly  enough  space  and  contain  their  integer 
number  of  rows  and  columns  as  the  first  two  data  values.  Subroutine  JACOBI  is  called 
on. 

Calling  Sequence— M0DES(A(1C),B(IC),Z(IC),R(IC),C(IC)) 

Subroutine  MASS 

Purpose— If  a  dynamic  vibration  problem  is  referred  to  a  set  of  coordinates  consist¬ 
ing  of  the  deflections  and  the  rotations  0i  at  N  collocation  points  along  the  beam 
under  consideration,  then  this  subroutine  generates  the  2N  by  2N  inertia  matrix  [A] 
which  appears  in  the  following  expression  for  kinetic  energy. 


Restrictions— The  mass  and  inertia  data  inputs  to  this  subroutine  are  to  be  supplied 
as  piecewise  con*:v.uous  slices;  however,  these  arrays  may  be  of  arbitrary  size  and  differ¬ 
ent  in  length  from  each  other.  The  number  of  collocation  points  N  which  determines  the 
ultimate  size,  2N  by  2N,  of  the  output  inertia  matrix,  is  also  chosen  arbitrarily. 

Calling  Sequence— MASS(X(IC),DMPL(IC),RIPL(IC),CM(IC),A(IC)) 

Here 

X  is  an  N*1  matrix  of  collocation  points  referred  to  an  arbitrary  origin. 

DMPL  is  an  NDM*4  matrix  of  distributed  mass  per  unit  length  slices,  in  which 
Col  1  is  the  location  of  the  rear  of  a  slice. 

Col  2  is  the  location  of  the  front  of  a  slice. 

Col  3  is  the  mass  value  at  the  roar  of  the  slice. 

Col  4  is  the  mass  value  at  the  front  of  the  slice. 

RIPL  is  an  NRI*4  matrix  of  distributed  rotary  inertia  per  unit  length  slices.  The 
columns  here  are  similar  to  DMPL. 

CM  is  an  NCM*4  matrix  of  concentrated  mass  items,  where 
Col  1  is  the  attach  point  location  for  each  item. 

Col  2  is  the  mass  at  this  location. 

Col  3  is  the  location  of  its  center  of  gravity. 

Col  4  is  the  amount  of  inertia  about  the  center  of  gravity. 

A  is  a  2N*2N  output  inertia  matrix. 

NOTE:  Since  this  applies  to  DMPL,  KIPL,  and  CM,  the  location  of  the  values  may  not 
go  beyond  the  limits  of  the  collocation  points  in  either  direction. 
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Subroutine  STIFF 

Purpose— If  a  dynamic  vibration  problem  is  referred  to  a  set  of  coordinates  consist¬ 
ing  of  the  deflections  and  the  rotations  0t  at  N  collocation  points  along  the  beam 
under  consideration,  then  this  subroutine  generates  the  2N  by  2N  stiffness  matrix  [K] 
which  appears  in  the  following  expression  for  the  strain  energy 


Restrictions— The  stiffness  and  shear  data  inputs  to  this  subroutine  are  to  be  sup¬ 
plied  as  piecewise  continuous  slices;  however,  these  arrays  may  be  of  arbitrary  size  and 
different  in  length  from  each  other.  The  numbej  of  collocation  points  N,  which  deter¬ 
mine  the  ultimate  size  (2N  by  2N)  of  the  output  stiffness  matrix,  is  also  chosen  arbitrarily. 

Calling  Sequence— STIFF(X(IC),EI(IC),GA(IC),K(IC)) 
where 

X  is  an  N  by  1  matrix  of  collocation  points  referred  to  an  arbitrary  origin. 

El  is  an  NEI  by  4  matrix  of  bending  stiffness  slices,  where 
Col  1  is  the  location  of  the  rear  of  a  slice. 

Col  2  is  the  location  of  the  front  of  a  slice. 

Col  3  is  the  stiffness  value  at  the  rear  of  a  slice. 

Col  4  is  the  stiffness  value  at  the  front  of  a  slice. 

GA  is  an  NGA  by  4  matrix  of  shear  stiffness  slices,  where  the  columns  here  are 
similar  to  those  for  the  El  distribution. 

K  is  the  output  stiffness  matrix  size  2N  by  2N. 

NOTE:  Since  this  applies  to  El  and  GA,  the  location  of  the  values  may  not  go  beyond 
the  limits  of  the  collocation  points  in  either  direction. 


Subroutine  LIST 

Purpose— This  subroutine  prints  out  the  elements  of  a  matrix  (A]  and  identifies  each 
by  its  row  and  column  number.  The  user  must  supply  an  alphanumeric  name  ALP  and 
integer  number  NUM  to  identify  the  matrix.  This  is  to  maintain  consistency  with  sub¬ 
routines  FILE  and  CALL. 

Restrictions— The  matrix  must  have  its  integer  number  of  rows  and  columns  as  the 
first  two  data  values. 


Calling  Sequence  —  LIST(  A( IC ) , A  LP,N UM ) 
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Subroutine  PUNCH 

Purpose—' This  subroutine  punches  out  a  matrix  [A],  size  n*m,  one  column  at  a 
time  in  any  desired  format.  The  argument  FOR.  must  reference  a  FORTRAN  format 
statement  that  has  been  entered  as  a  positive  array.  It  must  include  the  outer  paren¬ 
thesis  but  not  the  word  FORMAT.  The  argument  HEAD  must  be  a  single  BCD  word 
used  to  identify  the  matrix.  Each  column  is  designated  and  restarts  use  of  the  FORMAT 
statement. 

This  subroutine  requires  n+3  dynamic  storage  locations. 

Restrictions— The  matrix  ( A J  must  have  exactly  enough  space  and  contain  the  inte¬ 
ger  number  of  rows  and  columns  as  the  first  two  data  values.  Punched  cards  must  be 
asked  for  on  the  job  request  form. 

Calling  Sequence— PUNCH(A(IC), HE  AD, FOR(IC)) 


Matrix  Data  Storage  and  Retrieval 

The  ability  to  store  and  retrieve  matrixes  from  tape  is  easily  achieved  through  the 
use  of  the  FILE  and  CALL  subroutines.  Matrixes  are  identified  by  an  alphanumeric 
name,  integer  problem  number,  and  the  core  address  of  or  for  the  matrix.  The  CALL 
subroutine  searches  the  matrix  storage  tape  on  logical  16  and  brings  the  desired  matrix 
into  core.  The  FILE  subroutine  writes  a  matrix  onto  the  logical  30  tape.  Subroutine 
END.MOP  causes  all  matrixes  from  the  logical  30  tape  to  be  updated  onto  the  logical  16 
tape.  In  case  of  duplicate  matrixes  the  one  from  logical  30  replaces  the  one  on  logical 
16.  A  matrix  which  has  been  filed  cannot  be  called  until  an  ENDMOP  operation  has  been 
performed.  To  create  a  new  tape  the  user  merely  sets  control  constant  NOCOPY  nonzero 
and  has  a  scratch  tape  mounted  on  logical  16.  The  user  should  check  the  section  on  con¬ 
trol  cards  and  deck  setup  to  determine  control  card  requirements. 


Subroutine  CALL  or  FILE 

Purpose— To  allow  the  user  to  retrieve  or  store  matrixes  on  magnetic  tape,  see  above. 
The  H  argument  must  be  a  6-<:haracter  alphanumeric  word  and  N  must  be  an  integer 
number,  both  of  which  are  used  to  identify  the  matrix. 

Restrictions— See  above.  The  matrix  must  have  exactly  enough  space  and  contain 
the  integer  number  of  rows  and  columns  as  the  first  two  data  values. 

Calling  Sequence-C ALL  (H,N,A(IC))  or  FILE  (A(IC),H,N) 

Subroutine  ENDMOP  or  LSTAPE 

Purpose— Subroutine  ENDMOP  should  be  used  in  conjunction  with  subroutines  CALL 
and  FILE,  see  above.  It  causes  matrixes  which  have  been  filed  by  FILE  on  logical  30  to 
be  updated  onto  logical  16.  A  cail  to  subroutine  LSTAPE  will  cause  the  output  of  the 
name,  problem  number,  and  size  of  every  matrix  stored  on  tape  on  logical  16. 

Restrictions— See  above. 

Calling  Sequence  —  EN  DMOP  or  LSTAPE 
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Subroutine  SIMEQN 


Purpose— This  subroutine  solves  a  set  of  up  to  10  linear  simultaneous  equations  by 
the  factorized  inverse  method.  The  pi^blem  size  and  all  input  and  output  values  are 
communicated  as  a  single,  specially  formatted,  positive  input  array.  The  array  argument 
must  address  the  matrix  order  (N)  which  is  input  by  the  user.  The  first  data  value  must 
be  the  integer  order  of  the  set  (or  size  of  the  square  matrix)  followed  by  the  coefficient 
matrix  (A]  in  column  order,  the  boundary  vector  {  b}.  and  space  for  the  solution  vec¬ 
tor  {s}: 

[A]  {S}  -  {B}. 

Restrictions— The  integer  count  and  matrix  size  must  be  integers;  all  other  values 
must  be  floating  point.  The  coefficient  matrix  is  not  modified  by  SIMEQN.  Hence, 
changes  to  {b}  only  allow  additional  solutions  to  be  easily  obtained. 

Calling  Sequence— SIMEQN(A(N)) 


where  the  array  is  formatted  exactly  as  follows: 


IC,N,A(1,1),A(1,2), 


A(N,N),B1, . . . ,  BN,S1 . SN 


Subroutine  LSTSQU 
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Purpose— This  subroutine  performs  a  least  squares  curve  fit  to  an  arbitrary  number  f 

of  X,  Y  pairs  to  yield  a  polynomial  equation  of  up  to  order  10.  Rather  than  using  a  | 

double  precision  matrix  inverse,  this  subroutine  calls  on  subroutine  SIMEQN  to  obtain  a  | 

simultaneous  solution.  | 

I 
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This  subroutine  requires  2*M  dynamic  storage  core  locations. 

Restrictions— All  values  must  be  floating-point  numbers  except  N  and  M,  which  must 
be  integers.  N  is  the  order  of  the  polynomial  desired  and  is  one  less  than  the  number  of 
coefficients  desired.  M  is  the  array  length  of  the  independent  X  or  dependent  Y  values. 

Calling  Sequence— LSTSQU(N,M,X(DV),Y(DV),a(DV)) 


Subroutine  IRRADI  or  IRRADE 


Purpose— These  subroutines  simulate  a  radiosity  network*  within  a  multiple  gray 
surface  enclosure  containing  a  nonabsorbing  media.  The  input  is  identical  for  both  sub¬ 
routines.  However,  IRRADE  utilizes  explicit  equations  to  obtain  the  solution  by  relaxa¬ 
tion,  and  IRRADI  initially  performs  a  symmetric  matrix  algebra  inverse  and  thereafter 
obtains  the  exact  solution  implicitly  by  matrix  multiplication.  The  relaxation  criteria  of 
IRRADE  is  internally  calculated  and  severe  enough  so  that  both  routines  generally  yield 
identical  results.  However,  IRRADE  should  be  used  when  temperature-varying  emissivi- 
ties  are  to  be  considered,  and  IRRADI  should  be  used  when  the  surface  emissivities  are 
constant.  Both  subroutines  solve  for  the  J-node  radiosity,  obtain  the  net  radiant  heat 
flow  rates  to  each  surface,  and  return  them  sequentially  in  the  last  array  that  was  initially 
used  to  input  the  surface  temperatures.  The  user  need  not  specify  any  radiation  con¬ 
ductors  within  the  enclosure. 

Restrictions— The  Fahrenheit  system  is  required.  The  arbitrary  number  of  tempera¬ 
ture  arguments  may  be  constructed  by  a  preceding  BLDARY  call.  The  emissivity,  area, 
temperature-Q  and  upper  half-FA  arrays  must  be  in  corresponding  order  and  of  exact 
length.  The  first  data  value  of  the  FA  array  must  be  the  integer  number  of  surfaces  and 
..  j  second  the  Stephan-Boltzmann  constant  in  the  proper  units  and  then  the  FA  floating¬ 
point  values  in  row  order.  The  diagonal  elements  (even  if  zero)  must  be  included.  As 
many  radiosity  subroutine  calls  as  desired  may  be  used.  However,  each  call  must  have 
unique  array  arguments.  The  user  should  fellow  the  radiosity  routine  by  SCALE, 
BRKARY,  or  BKARAD  to  distribute  the  Q’s  to  the  proper  source  locations. 

Calling  Sequences— IRRADI(AA(IC),Ae(IC),AFA(IC),ATQ(IC))  or 
IRRADE(AA(IC),Ae(IC),AFA(IC),ATQ(IC)) 

The  arrays  are  formatted  as  follows: 


AA(IC),AI,A2,A3,A4, ....  AN, END 
Ae(IC),el,e2,e3,e4, ....  eN.END 

AFA(IC),N,o,FA(l,l  ),FA(1,2),F A(1,3),F A(1,4),FA(1,5), . . .  FA(1,N) 
FA(2,2),FA(2,3),FA(2,4  ),FA(2,5), . . .  FA(2,N) 

FA(N-2,N-2),KA(N~2,N-1),FA(N-2,N) 

FA(N-1,N-1),FA(N-1,N) 

FA(N  N),END 

ATQ(IC),T1,T2,T3, . . .  TN,END 


*A.  K.  Oppcnhcim,  "Radiation  Analysis  by  the  Network  Method,”  Trans.  ASME  78,  725-735  (1956). 
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where  FA(1,2)  is  defined  as  A(l)  *  F(l,2).  After  the  subroutine  has  been  performed  the 
ATQ  array  is  ATQ(1C),Q1,Q2,Q3, . . .  QN.END. 

Since  FAj(l,2)  =  FA2(2,1)  only  the  upper  half  triangle  of  the  full  FA  matrix  is 
required.  1RRAD1  inverts  this  half-matrix  in  its  own  area;  hence,  approximately  300 
surfaces  may  be  considered  using  C1NDA  on  a  65k-core  machine. 


Subroutine  SLRAOI  or  SLRADE 

Purpose— These  subroutines  are  very  similar  to  IRRAD1  and  IRRADE  but  are  de¬ 
signed  to  solve  for  the  solar  heating  rates  within  an  enclosure.  SLRADI  inverts  half  a 
symmetric  matrix  in  order  to  obtain  implicit  solutions  while  SLRADE  obtains  solutions 
explicitly  by  relaxation.  SLRADE  should  be  used  when  temperature  varying  solar  emissivi- 
ties  are  to  be  considered.  The  second  data  value  of  the  AFA  array  must  be  the  solar  con¬ 
stant  in  the  proper  units.  The  AT  array  allows  the  user  to  input  the  angle  (degrees)  be¬ 
tween  the  surface  normal  and  the  surface-sun  line.  The  AI  array  allows  the  user  to  input 
an  illumination  factor  for  each  surface  which  is  the  ratio  from  zero  to  one  of  the  unshaded 
portion  of  the  surface.  The  solar  constant  S,  AT,  and  AI  values  may  vary  during  the 
transient  for  both  routines.  No  input  surface  temperatures  are  required.  The  absorbed 
heating  rates  are  returned  sequentially  in  the  AQ  array;  the  user  may  utilize  SCALE, 
BRKARY,  or  BKARAD  to  distribute  the  heating  rates  to  the  proper  source  locations. 

Restrictions— These  routines  are  independent  of  the  temperature  system  being  used. 
All  of  the  array  arguments  must  reference  the  integer  count  set  by  the  CINDA  preproc¬ 
essor  and  be  of  the  exact  required  length.  As  many  calls  as  desired  may  be  made,  but 
each  call  must  have  unique  array  arguments. 

Calling Se<7uences-SLRADI(AA(IC),Ae(IC),AFA(IC),AT(IC),AI(IC),AQ(IC))  or 
SLRADE(AA(IC),Ae(IC),AFA(IC),AT(lC),AI(IC),AQ(IC)) 


Subroutine  SCRPFA 

Purpose— To  obtain  the  script  FA  value  for  radiant  transfer  within  an  enclosure. 

The  input  arrays  are  formatted  as  shown  for  subroutines  IRRADI  and  IRRADE.  The 
second  data  value  in  the  AFA  array  is  used  as  a  final  multiplier.  If  1.0  the  script  FA 
values  are  returned;  if  o  then  script  o  FA  values  are  returned.  The  script  FA  values  are 
returned  in  the  ASFA  ar.-ay  which  is  formatted  identically  to  the  AFA  array  and  may 
overlay  it. 

Restrictions— All  array  arguments  must  reference  the  integer  count  set  by  the 
CINDA  preprocessor,  and  all  arrays  must  be  exactly  the  required  length. 

Ca//mgSequence— SCRPFA(  A  A(IC),Ae(IC),AFA(IC),ASFA(lC)) 

NOTE:  Subroutine  SYMLST(ASFA(1C)+3,ASFA(IC)+1)  may  be  called  to  list  the  matrix 
values  and  identify  them  by  row  and  column  number.  This  routine  and  the  implicit 
radiosity  routines  finalize  the  half-symmetric-coefficient  matrix  and  call  on 
SYMINV(AFA(IC)+3.AFA(IC)+1)  to  obtain  the  symmetric  inverse. 
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Subroutine  ABLATS 

Purpose— ABLATS  provides  a  simple  ablation  (sublimation)  capability  for  the  CINDA 
user.  The  user  constructs  the  three-dimensional  network  without  considering  the  ablative. 
Then  in  Variables  2  he  simulates  one-dimensional  ablative  attachments  by  calling  aBLATS. 
ABLATS  constructs  tne  one-dimensional  network  and  solves  it  by  implicit  forward- 
backward  differencing  (Crank-Nichoison  method)  using  the  time  step  set  by  the  execu¬ 
tion  subroutine.  Separate  ablation  arrays  (AA)  must  be  used  for  each  ABLATS  call.  Re¬ 
quired  working  space  is  obtained  Lorn  unused  program  COMMON.  Several  ABLATS 
calls  thereby  share  unused  COMMON.  The  user  must  call  subroutine  PNTABL(AA)  in  the 
OUTPUT  CALLS  to  obtain  the  ablation  totals  and  temperature  distribution. 

Restrictions— ABLATS  must  be  called  in  Variables  2  and  may  be  used  with  any 
execution  subroutine.  Subroutines  D1DEG1,  NEWTR4,  and  INTRFC  are  called.  All 
units  must  be  consistent.  The  Fahrenheit  system  is  required,  femperature-varying  mate¬ 
rial  property  arrays  must  not  exceed  60  doublets.  Bivariate  material  properties  may  be 
simulated  by  calling  BVSPSA  prior  to  ABLATS.  Cross-sectional  area  is  always  considered 
unity.  Thermal  conductivity,  Stephan-Boltzmann  constant,  and  density  units  must  agree 
in  area  and  length  units. 

This  subroutine  requires  3*(NSI  vl)  d:  namic  storage  core  locations. 

Calling  Sequence- ABLATS(AA{\C),R,CP,G,T,C) 
where 

C  is  the  capacitance  location  of  the  three-dimensional  node. 

T  is  the  temperature  location  of  the  three-dimensional  node. 

G  is  the  location  of  the  material  thermal  conductivity  or  the  starting  location  (inte¬ 
ger  couni ,  of  a  doublet  G  vs  T  array. 

CP  is  the  location  of  the  material  specific  heat  or  the  starting  location  (integer 
count)  of  a  doublet  Cp  vs  T  array. 

R  is  the  location  of  the  material  density  or  the  starting  location  (integer  count)  of  a 
doublet  p  vs  T  array. 

AA(IC)  is  the  starting  location  of  the  ablation  array  which  must  be  formatted  as 
follows: 

AA(IC)+1  -the  ablative  link  number,  a  user-specified  identification  integer. 

AA(lC)+2--integer  number  of  sublayers  (NSL)  desired;  ABLATS  subtracts 
from  this  the  number  of  sublayers  ablated. 

AA(lC)+3— the  initial  temperature  of  the  material;  ABLATS  replaces  this  with 
the  outer  surface  temperature,  always  in  degrees  F. 

AA(IC)+4— the  impressed  outer  surface  heating  rate  per  unit  area,  radiation 
rates  not  included. 

A A(IC)+5— material  thickness;  this  is  replaced  by  the  sublayer  thickness. 

AA(IC)+6— surface  area  of  the  three  -dimensional  node;  need  not  be  unity. 

AA(IC)+7— ablation  temperature,  degrees  F. 

AA(IC)+8— heat  of  ablation. 

AA(lC)+9— Stephan-Boltzmann  constant  in  consistent  units. 

AA(1C)+10— surface  emissivity. 
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AA(IC)+11— space  “sink”  temperature,  degrees  F. 

AA(IC)+12-SPACE,N,END  where  N  equals  NSI  +  4. 

NOTE:  The  outer  surface  radiation  loss  is  integrated  over  the  time  step. 

Subroutine  LQDVAP 

Purpose— This  subroutine  allows  the  user  to  simulate  the  addition  of  liquid  to  a 
node.  The  network  data  is  prepared  as  though  no  liquid  exists  at  the  node  and  is  solved 
that  way  by  the  network  execution  subroutine.  Then  LQDVAP,  which  must  be  called 
from  Variables  2,  corrects  the  nodal  solution  in  order  to  account  for  the  liquid.  If  the 
nodal  temperature  exceeds  the  boiling  point  of  the  liquid,  it  is  set  to  the  boiling  point. 

The  excess  energy  above  that  required  to  reach  the  boiling  point  is  calculated  and 
considered  as  absorbed  through  vaporization.  If  the  liquid  is  completely  vaporized  the 
subroutine  deletes  its  operations.  The  method  of  solution  holds  very  well  for  explicit 
solutions,  but  may  introduce  some  error  when  large  time  steps  are  used  with  implicit 
solutions. 

Restrictions— This  subroutine  must  be  called  from  Variables  2. 

Calling  Sequence— LQDVAP(T,C,A(  IC) ) 
where 

T  is  the  temperature  location  of  the  node. 

C  is  the  capacitance  location  of  the  node. 

A  +  1  contains  the  initial  liquid  weight. 

A  +  2  contains  the  liquid  specific  heat. 

A  +  3  contains  the  liquid  vaporization  temperature. 

A  +  4  contains  the  liquid  heat  of  vaporization. 

A  +  5  receives  the  liquid  vaporization  rate  (weight/lime). 

A  +  6  receives  the  liquid  vaporization  total  (total  weight). 

A  +  7  contains  the  liquid  initial  temperature. 

Subroutine  BIV  LV 

Purpose— This  subroutine  allows  the  user  to  specify  the  percentage  flow  rates 
through  two  parallel  tubes  with  common  end  points.  One  tube  must  consist  of  a  single 
flow  conductor  (Gl)  while  the  other  tube  may  consist  of  one  or  more  sequential  flow 
conductors  (G2(I),  I  =  1,N).  The  ratio  of  flow  through  Gl  divided  by  the  total  flow 
may  be  calculated  ir.  any  desired  manner  and  must  be  supplied  as  the  argument  W.  The 
conductor  values  of  either  one  tube  or  the  other  are  reduced  in  order  to  achieve  the  de¬ 
sired  percentage  flow  rates  irregardless  of  the  pressure  drop. 

Restriction:  -N  must  be  an  integer.  G2  must  address  the  first  of  the  sequential 
conductors  in  that  tube. 

Calling  Seq ucncc —BIV LV ( N.W.G  1  ,G2( D V ) ) 
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Subroutine  LINE 

Purpose— This  subroutine  computes  the  steady  state  changes  ir  the  thermodynamic 
and  flow  properties  through  a  line  of  length  L.  The  upstream  properties  must  be  de¬ 
fined  and  supplied.  The  following  equations  are  simultaneously  solved  in  an  iterative 
fashion: 

puvu  =  pdvd  (one-dimensional  conservation  of  mass), 
where  u  denotes  upstream  and  d  denotes  downstream; 

,  ,  v“2  .  Q  .  *  vd2  .  . 

lu  +  ~n~  +  ttt  =  L  +  ~x~  (one-dimenstonai  energy 

^  w  i  ..  . 

equation); 

Pu  ~  Pd  (momentum  equation,  simplified  because  of  a 
very  small  pressure  drop  and  low  velocities); 

P  =  f(X,P,h) 

T  =  g(X,P,h)  (equations  of  state); 

XL  -  h(X,P,h) 

Q  =•  h(Wp)L(Tw  -  T)  (energy  loss) 

Re  =  W(4A/Wp)/Ap  (Re  is  the  Reynolds  number) 

Pr  =  Cp  p/ii  (Pr  is  the  Prandtl  number) 
h  =  0.332  (Re°-5)(Pr0-333)  for  Re  <  2100  (laminar  flow) 
h  =  0.023  (Re°-8)(Pr0-4)  for  Re  >  2100  (turbulent  flow) 

Cp,  p,  and  K  are  obtained  from  subroutine  TRNPRT. 

Restrictions— Subroutine  LINE  assumes  that  the  flowing  fluid  is  composed  of  a  per¬ 
fect  noncondensable  gas  and  a  perfect  condensable  gas.  This  assumption  involves  the 
STATE  subroutine  which  ia  called  on.  However,  LINE  does  not  need  the  variables  XL 
and  T  to  evaluate  the  transport  properties  and  the  boat  transfer  coefficient  h  for  the  cal¬ 
culation  of  Q.  The  thermodynamic  property  arguments  are  upstream  properties  when 
calling  LINE,  but  the  downstream  thermodynamic  properties  are  in  the  same  locations. 

Calling  Sequence— 

LINE(A,Wp,L,Tw,W,Al(IC),A2(IC),A3(IC), . . .  ,A10(IC),P,T,XL,X,I,V,p,Q) 

where 

.a  -  flow  are? 

Wp  =  wetted  perimeter 
L  =  tube  length  (floating  point) 

Tw  =  wall  temperature 
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VV  =  mass  flow  rate 

A1  =  the  doublet  interpolation  array  of  condensable  gas,  p  vs  T 
A2  =  the  doublet  interpolation  array  of  noncondensable  gas,  p  vs  T 

A3  =  the  doublet  interpolation  array  of  condensable  gas,  k  vs  T 

A4  =  the  doublet  interpolation  array  of  noncondensable  gas,  k  vs  T 

A5  =  the  doublet  interpolation  array  of  condensable  gas,  Cp  vs  T 

A6  =  the  doublet  interpolation  array  of  noncondensablc  gas,  Cp  vs  T 

A10  =  the  doublet  interpolation  array  of  condensable  gas,  heat  of  vaporization  vs  T 
P  -  pressure 
T  =  temperature 
Xl  =  mass  fraction  of  liquid 
X  =  mass  fraction  of  noncondensable  gas 
I  =  enthalpy  (floating  point) 

V  =  v  =  velocity 
p  =  density 

Q  =  energy  lost  to  wall 


Subroutine  STATE 

Purpose  -This  subroutine  computes  the  thermodynamic  state  for  a  mixture  of  an 
assumed  noncondensable  gas  (hydrogen)  and  a  condensable  gas  v water  vapor).  The  sub¬ 
routine  establishes  whether  the  mixture  is  superheated  or  saturated  and  gives  its  density 
(pm),  temperature  (T),  and  liquid  mass  fraction  (XL).  The  hydrogen  mass  fraction  (Xh), 
mixture  pressure  (Pm),  and  mixture  enthalpy  (Im)  are  input.  Vapor  components  are 
assumed  to  be  perfect  gases;  that  is. 


Iv  -  CpvT  Ih  =  CphT 

where  subscripts  h  and  v  refer  to  hydrogen  and  water  vapor,  respectively.  The  liquid 
constituent  is  assumed  to  have  the  following  properties: 

oL  =  62.4  and  IL  =  Iv  -  HV, 

where 

HV  =  heat  of  vaporization 
Cp  -  specific  heat  at  constant  pressure 
R  =  gas  constant. 

If  the  mixture  is  saturated,  Pv  is  related  to  T  by  the  saturation  equation  of  sub¬ 
routines  PSOFTS  and/or  TSOFP.  Mixture  properties  are  obtained  from  the  following 
equations: 
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_ IX) _ 

Pm  =  Xh  1  (1  —  xh  —  XL)  |  XL 

Ph  Py  Pl 

and 

Im  "  XhIh+(l-XL)Iv-XL(HV). 

Restrictions— The  restrictions  are  those  that  are  imposed  for  perfect  gases  and  incom¬ 
pressible  liquids.  The  pressures  must  be  well  below  the  critical  point.  A  is  a  doublet  in¬ 
terpolation  array  of  HV  as  a  function  of  temperature  T.  Im  must  be  a  floating-point 
number. 

Calling  Sequence— STATE(X,P,I,p,XL,T,A(IC)) 


Subroutine  PSOFTS 

Purpose— This  subroutine  computes  the  saturation  pressure  of  water  vapor  as  a 
function  of  gas  temperature.  The  relationship  used  is 

Pc  x  (a  +  Bx  +  CxA 

loe>0  7  ■  tX  TTdx  r 

where  x  =  Tc  -  Ts,  A,  B,  C,  and  D  are  constants,  and  Pc  and  Tc  are  critical  points. 
Restrictions— The  gas  temperature  should  be  between  10°  and  150°  C. 

Calling  Sequence- PSOFTS(TS.P) 


Subroutine  TSOFP 

Purpose —This  subroutine  computes  gas  temperature  as  a  function  of  the  saturation 
pressure  of  water  vapor,  using  the  same  relation  as  in  PSOFTS. 

Restrictions— The  saint  restrictions  apply  to  TSOFP  as  to  PSOFTS. 

Calling  Sequence— TSOFP(P,TS) 


Subroutine  TRNPRT 

Purpose— This  subroutine  calculates  the  transport  properties  of  a  two-component 
gas  mixture. 

Restrictions— Only  a  two-component  gas  mixture  is  allowed,  and  the  component 
properties  must  have  already  been  evaluated  at  the  desired  temperature. 
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Calling  Sequence  — TRNPRT(  VI  ,V2,G  1  ,G2,C1  ,C2,V ,G ,C,P1 ) 
where 

VN  is  the  viscosity  of  component  N. 

GN  is  the  thermal  conductivity  of  component  N. 

CN  is  the  specific  heat  of  component  N. 

PI  is  the  percent  (by  weight)  of  component  1. 

V,  G,  and  C  are  the  viscosity,  thermal  conductivity,  and  specific  heat  of  the  mixture. 


Subroutine  CSGDMP  or  LSTPCS  or  QMAP 

Purpose— These  routines  are  designed  to  aid  in  the  checkout  of  thermal  problem 
data  decks  by  listing  the  pseudo-compute  sequence.  CSGDMP  calls  upon  Variables  1  and 
then  prints  out  each  relative  diffusion  node  number  with  the  capacitance  and  CSGMIN 
value  of  the  node.  For  each  node,  all  three  routines  identify  the  attached  conductors  by 
relative  conductor  number  and  type,  and  by  the  relative  number  of  the  adjoining  node. 
CSGDMP  also  lists  the  conductance  of  the  attached  conductor  and  the  type  of  the  ad¬ 
joining  node.  Either  the  SPCS  or  the  LPCS  option  may  bs  used.  While  the  LPCS  option 
allows  every’  conducior  attached  to  a  node  to  be  identified,  the  SPCS  option  identifies 
only  conductors  for  the  first  relative  node  number  on  which  they  occur.  After  the  dif¬ 
fusion  nodes  are  processed,  the  connection  information  for  the  arithmetic  nodes  is 
listed.  After  listing  the  above  information,  control  passes  to  the  next  sequentially  listed 
subroutine. 

QMAP  has  all  the  properties  of  CSGDMP.  In  addition,  it  prints  the  temperatures  of 
each  node  and  adjoining  node,  and  the  flux  between  them. 

Restrictions— These  routines  are  generally  called  from  EXECTN.  CSGDMP  and 
QMAP  should  never  be  called  from  Variables  1. 

Calling  Sequence— CSGDMP  or  LSTPCS  or  QMAP 


Subroutine  TSAVE 

Purpose— This  subroutine  ger;  'rates  an  external  plotting  data  output  file  (unit  24) 
that  can  be  used  with  the  external  plowing  option  to  plot  nodal  temperature  vs  time. 
TSAVE  records  each  nodal  temperature  at  TIMEN  and  also  saves  the  actual  node 
numbers. 

Restrictions— This  routine  should  not  be  called  more  than  2000  times. 


Calling  Sequence— TSAVE 
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Internal  Subroutines 


Name 

Name 

Name 

BIT 

ITRATE 

SMOPAS 

COPY 

POLYADD 

TOPLIN 

GENM 

PYMLT1 

UNPAK 

GENST 

RTPOLY 

UPDMOP 

HEDCOL 

SETUP 

WRTARY 

1NPUTG 

INPUTT 

SKPLIN 

WRTL08 

These  subroutines  are  called  from  other  routines  and  not  normally  by  the  user. 
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Appendix  A 

SAMPLE  PROBLEM  1A 


ORIGINAL  RUN 

A  perfectly  insulated  one-dimensional  bar  has  a  constant  heating  rate  applied  to  one 
end.  Obtain  the  10-min  transient  temperature  response,  at  half-minute  intervals,  of  the 
bar  ends  and  at  points  1/4,  1/2,  and  3/4  of  the  way  along  the  bar.  The  bar  is  initially 
at  80°F  and  receives  a  constant  heating  rate  of  3.0  Btu/min.  The  length  of  the  bar  is 
4  in.,  and  it  has  a  cross-sectional  area  of  1  sq.  in.  It  has  the  following  material  properties: 

density  =  172.8  lb/ft3 

specific  heat  =  0.35  Btu/lb°F 

thermal  conductivity  =  0.2  Btu/in.-min.-°F 

Figure  Ala  shows  a  schematic  of  the  physical  problem  with  the  nodes  appropriately 
placed  and  the  dashed  lines  indicating  the  lumping  of  the  system  for  capacitance  purposes. 
The  network  representation  is  illustrated  in  F;^.  Alb. 


(a) 


(a)  Parameter  lumping 


T1  T2  T3  T4  T5 


(b)  One-dimensional  network 


Fig.  A1  —Representation  of  a  perfectly  insulated  bar 

Capacitor?  receive  the  same  number  as  the  temperatures  but  with  a  C  prefix.  From 
the  above  information,  we  immediately  calculate 

C2  =  C3  =  C4  =  p*V*Cp  =  0.035  Btu/°F 

Cl  =  C5  =  C2/2.0  =  0.0175  Btu/°F 

G1  =  G2  =  G3  =  G4  =  k*Ac/C  =  0.2  Btu/min°F, 
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where  V  =  C*Ac;  length  times  cross-sectional  area. 


Since  this  is  not  a  RECALL  run,  the  first  data  card  should  be  blank. 


To  apply  explicit  forward  differencing  to  this  problem,  we  must  utilize  the  CNFRWD 
execution  subroutine  which  requires  the  short  pseudo-compute  sequence.  Hence,  the  title 
block  is  as  follows: 

8 

1 

BCD  3 THERMAL  SPCS 

BCD  9 SAMPLE  PROBLEM  NO. IA 

END 


The  nodal  block  is  next  and  requires  the  node  number,  initial  temperature,  and  capaci¬ 
tance  of  each  node  be  listed. 

8 

i 

BCD  3N0DL  DATA 

I ,80., .0175, 5, 80., .01 75 
GEN  2,3,1,80.,. 035 , I . , t  * , 1 . 

END 


The  conductor  block  requires  that  each  conductor  number  be  listed  with  the  node  num¬ 
bers  at  either  end,  and  the  conductor  value. 

8 

i 

BCD  3CQN0UCT0R  DATA 

GEN  1,4, I, 1,1, 2,1, .2,1. ,!.,!. 

END 

The  only  control  constants  required  for  CNFRWD  are  as  follows: 

8 

i 

BCD  3C0NSTANTS  DATA 

TIMEND, 10. , OUTPUT , . 5, CSGFAC.2 . 

END 

There  are  no  array  data  and  only  one  execution  call;  hence, 
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18  21  25 

11  11 
BCD  3 ARRAY  DATA 
bND 

F  DIMENSION  X(  100) 

BCD  3EXECUTI0N 
F  NDIM=1 00 

F  NTH=0 

CSUDMP 

CNFR1VD 

END 

There  are  no  second  variables  operations,  but  we  must  apply  the  heating  rate  in  the  first 
variables; 


BCD  3VARI  ABLliS  I 
STFSbP  (3. ,01 ) 

bND 

BCD  3VARI ABLES  2 
END 

The  following  completes  the  data  input. 


BCD  3GUTPUT  CALLS 
PRNTiMP 

END 

BCD  3END  UF  DATA 

Since  PRNTMP  lists  the  relative  node  numbers,  and  not  the  actual  ones,  the  node  dic¬ 
tionary  will  have  to  be  consulted  for  conversion  of  relative  to  actual. 

The  above  problem  data  deck  processed  by  the  C1NDA  program  on  the  CDC-3800 
as  a  standard  run  produces  the  output  as  given  in  the  following  printouts. 

NOTE:  The  only  alternative  to  the  BCD  3END  OF  DATA  card  is  a  parameter  change. 
A  new  job  would  require  another  set  of  control  cards. 


SEQUENCE  01713  STARTEO  PRINTING  07/13/73  A!  185705  ON  LP01 

DRUM  SCOPE  2.1  COMPUTER  TWO.  MAX.  OEMANO  IS  57000B  VERSION  00A  10/31/72 

SEQUENCE  NUMBER  001713  STARTED  AT  TIME  1B2721  OATEO  07/13/73 
JOB (5) • *0800800. 89RRCC. 5 

EQUIP. 09«ICIN0A  MASTER. 1.1. 999) .RO.Hj.OA 
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SAMPLE  PROBLEM  NO.  1A  (Continued) 

TSAVE  AND  PLOT  RUN 

This  is  an  example  of  a  TSAVE  run  that  was  made  after  the  problem  had  been  satis¬ 
factorily  debugged  and  run.  Since  that  run  had  produced  printed  output,  all  calls  to  out¬ 
put  subroutines  were  removed  from  the  deck  to  save  processing  time.  Hence,  the  Execu¬ 
tion  block  is  now 


1 

7 

21 

25 

1 

1 

4 

4 

F 

DIMENSION 

X( 

100) 

BCD  3EXECUTI0N 


F  NDIM  =  100 

F  NTH  =  0 

END 

TSAVE  is  the  only  subroutine  in  Output  Calls. 

8  12 

1  I 

BCD  30UTPUT  CALLS 
TSAVE 

bND 

Since  the  time  and  temperature  limits  were  determined  from  the  output  in  the  previ¬ 
ous  run,  it  is  possible  to  supply  the  plotting  data  so  that  the  plot  program  can  be  run 
immediately  after  the  C1NDA  problem  (in  the  same  job).  This  eliminates  the  need  for 
equipping  the  TSAVE  output  tape  (tape  24).  The  plotting  Hata  are  as  follows: 


1 

12 

22 

32 

42 

4 

4 

4 

i 

'T 

4 

CINDA 

SAMPLE 

0.00 

PROBLEM  1A 
10.00 

0.00 

320.00 

card  3  (blank) 

HOF 

The  above  plot  data  and  revised  problem  data  deck  for  problem  1 A  produce  the 
following  output  and  plots  when  processed  on  the  CDC-3800  and  plotted  on  the 
CalComp  5G5  plotter  (the  actual  dimensions  of  the  X  and  Y  axes  are  7  and  9  in.,  re¬ 
spectively).  Plots  of  the  data  (Fig.  A2a-e)  follow  the  printouts. 
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temperature  profiles  for  CINDA  Sample  Problem  1A.  Plots  are  shown  in  the  order  in  which  they  are  plotted.  Time  and 
temperature  numbers  will  be  in  the  B  format,  rather  than  in  the  F  format  as  shown. 
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Appendix  B 

SAMPLE  PROBLEM  IB 


ORIGINAL  RUN 

Sample  problem  1A  was  linear  and  can  be  rigorously  solved  by  means  of  the  Laplace 
transform.  However,  the  introduction  of  nonlinearities  makes  rigorous  solutions  virtually 
impossible  and  makes  the  use  of  finite  difference  techniques  mandatory.  To  demonstrate, 
apply  the  following  nonlinearities  to  sample  problem  1A  and  obtain  the  solution. 

1.  Both  ends  of  the  bar  are  uninsulated  and  allowed  to  radiate  to  absolute  zero. 

The  Stephan-Boltzmann  constant  is  a  =  1.991E-13  Btu/min-in.2R°4,  and  the  emissivity 
varies  linearly  with  temperature  as  follows: 

e  =  0.4  at  -100°F 
e  =  0.8  at  300°  F. 

2.  The  thermal  conductivity  of  the  bar  varies  with  temperature  as  follows: 

k  =  0.15  at  -100°F  (Btu/in.-min-°F) 
k  =  0.25  at  100°F 
k  =  0.40  at  200° F 
k  =  0.60  at  300° F. 

3.  The  density  remains  unchanged  but  the  specific  heat  varies  with  temperature 
as  follows: 


Cp  =  0.3  at  -100° F  (Btu/lb-°F) 

=  0.39  at  100°F 
=  0.49  at  200°F 
=  0.65  at  300°F. 

4.  The  heating  rate  is  a  function  of  time  as  follows: 

q  =  3.0  at  0  min  (Btu/min) 
q  =  4.0  at  3  min 
q  =  4.0  at  7  min 
q  =  3.0  at  10  min. 

In  addition,  obtain  the  rate  of  heat  loss  and  the  integral  of  the  radiation  transfer 
from  the  unheated  end  of  the  bar.  The  network  representation  of  this  problem  (shown 
in  Fig.  Bl)  differs  only  slightly  from  problem  1A.  Now  however,  the  capacitances  are  a 
function  of  temperature.  We  therefore  require  multiplying  factors  such  that 
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Fig.  Bl  —Network  of  a  nonlinear  bar 


C  =  p  VCp(T),  MF  =  p  V 
MF  =  0.1  for  capacitors  2,  3.  and  4 
MF  =  0.05  for  capacitors  1  and  5. 

The  conductors  are  now 

G  =  k(Tm)Ac/t,  MF  =  Ac/H,  where  Tm  is  the  mean  of  the 
end  T’s. 

MF  =1.0  for  conductors  1,  2,  3,  and  4. 


A  radiation  conductor  requires  the  input  value  oeFA;  however,  FA  =  1.0,  hence 

Grad  =  oe(T| 

MF  =  1.991  Btu/min°F. 


Also, 


q  =  q(T|. 

The  capacitors  and  conductors  will  be  specified  with  CGS  and  CGD  calls. 


A.  Original  Run 

Since  this  is  not  a  RECALL  problem,  the  first  card  of  the  problem  data  deck  will 
be  blank.  The  rest  of  the  deck  may  be  constructed  as  follows: 
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8  12 
1  i 

BCD  3THERMAL  SPCS 
BCD  9 SAMPLE  PROBLEM  IB 
END 

BCD  3NODE  DATA 

CGS  1 ,8C. , A3, .05,2,80., A3,. 1,3, 80., A3,. I 
CGS  4,80. , A3, . I ,5,80. , A3, .05 
-10,-460. ,0 

END 

BCD  3 CONDUCTOR  DATA 

CGS  I , I ,2, A2, 1 . ,2,2,3, A2, ! . ,3,3,4, A2, 1 . ,4,4,b,A2, 1 . 

CGS  -1 1 , 1 ,!0,AI ,-l .99IE-! 3,-1 2,5, 1 0 ,A I ,-t .99IE-I3 
END 

BCD  3CUNSTANTS  DATA 

T I  MEND, 10. , OUTPUT, . 5, CSGFAC ,2 . , 4,0, 5, 0,6, STOR I ,7,STCR? 

BCD  3ARRAY  DATA 

1 , -100. ,.4, 300. ,0.8, END  $  EPSILON  VS  T 

2,  -100. ,.15,100. ,.25, 200. ,.4, 300., .6, END  $  K  VS  T 

3, -100. ,.3, 100. ,.39,200., .49,300. ,.65, END  $  CP  VST 

4,0. ,3. ,3. ,4. ,7. ,4. ,10. ,3., END  $  0  VS  TIME 
-5,QRATE,QT0TAL,END  $A  LABEL  ARRAY 

END 

F  DIMENSION  X(  100) 

BCD  3EXECUTI0N 
F  NDIM  =  100 

F  NTH  =  0 

STOREP  ( K6) 

CNFRWD 

F  IDCNT  =1  DC NT  +  1 

STURlP( K7 ) 

END 

BCD  3VARI ABLES  1 

D 1  DEG  I (TI ME M , A4 , 0 1 )  $  APPLY  HtATING  RATE 

END 

BCD  3 VAR I ABLES  2 

RDTNQSCTl  0,T5,G I  2 ,K4)  $UBTAI N  HEAT  FLOW  RATu 
0 INTcG( K4 , DTI MEU , K5)  $  INTEGRATE  SAME 

END 

BCD  3 OUTPUT  CALLS 
T PR  I  NT 

PPINTL  (A5,K4,K5) 

END 

BCD  3 END  OF  DATA 
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This  problem  will  be  stored  twice  on  tape  22.  The  original  data  will  be  stored  under 
the  I.D.  name,  ST0R1,  and  the  number,  0  (IDCNT).  The  final  values  will  be  identified 
as  ST0R2,  1  because  IDCNT  was  incremented.  (It  actually  would  not  have  been  neces¬ 
sary  to  use  IDCNT  in  this  case,  since  neither  call  to  ST0REP  was  in  a  loop.  The  second 
call  could  have  been  uniquely  identified  as  ST0R2,O.)  The  binary  constructed  subrou¬ 
tines  (processor)  will  be  stored  on  tape  40.  See  Section  VII  for  the  proper  deck  setup 
and  job  request  form. 

The  above  problem  data  deck  processed  by  the  CDC-3800  version  of  C1NDA  pro¬ 
duces  output  given  in  the  following  pages  (original  run). 
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SEQUENCE  NUMRER  004645  TERMINATED  AT  19204*  ELAP5E0  TIME*  00  MRS*  04  MIN*  31  SEC*  OR  271432  MILLISECONDS 


MARY  E.  GEALY 


The  stored  original  problem  of  sample  IB  is  used  to  illustrate  the  RECALL  option. 
The  initial  data  will  be  taken  and  given  a  new  TIMEND  of  5.0  min,  as  well  as  a  negative 
heating  rate. 

The  first  card  in  the  problem  data  deck  is 

1  13  22 

l  l  4 

RECALL  STORI  0 

The  title  block  is  changed  to 

8 

1 

BCD  3INTIAL  PARAMETERS 

BCD  3SAMPLE  PROBLEM  IB - RECALL 

END 

Since  there  are  no  temperature  or  conductor  data  changes,  the  nodal  and  conductor 
blocks  are  blank. 

8 

1 

BCD  3N0DE  DATA 
END 

BCD  3CUNUUCT0R  DATA 
END 

The  change  in  TIMEND  produces  the  following  constants  block: 

8 

l 

BCD  3 CONSTANTS  DATA 
TIMEND, 5.0 

END 

Array  4  is  changed  to  induce  a  negative  heating  rate. 

8 

I 

BCD  3ARRAY  DATA 

4,0., -3. ,3. ,-4 . ,7 . < -4. , 10., -3., END  $-0  VS  TIME 

END 

Since  no  operations  block  changes  are  allowed  in  parameter  runs,  the  data  deck  is  ter¬ 
minated  by 

8 

I 

BCD  3END  OF  DATA 

182 
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The  binary  program  tape  (processor)  was  stored  during  the  original  run,  so  recom¬ 
pilation  of  the  constructed  subroutines  is  not  necessary.  Since  the  STORE?  calls  are  still 
in  the  processor,  the  original  and  final  data  of  this  RECALL  problem  are  stored  on  drum 
unit  22.  I'he  unit  was  not  equipped  in  this  example,  however.  See  Section  Vll  for  the 
correct  deck  setup. 

The  processed  problem  printouts  of  the  second  run  follow. 
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1  5  12  22 
ill  i 

CINDA  SAMPLE  PROBLEM  IB 

0.00  10.00 

*  1 0* 

CINDA  SAME  PROBLEM  IB 

0.00  10.00 

+  10* 

EOF 


32 

1 

80.00 


42 

i 

?80.  00 


-470.00  -450.00 


The  printed  output  from  the  plot  run  are  as  follows. 


CINDA  SAMPLE  PROBLEM  18 

X*AXIS  LIMITS  ••  0.00*000f  1.00*091* 


T-AXIS  LIMITS  —  8. 00*001*  2J80*002 


all  nodes  hot  the  following  will  be  plotted  — 
10* 


MARY  E.  GEALY 
SAMPLE  NO.  IB  (Continued) 


TSAVE  AND  PLOT  RUN 

This  run  and  the  plots  are  similar  to  the  TSAVE  run  of  sample  1A.  The  following 
data  were  used. 


il 


i 


cinoa  sample  problem  ib 

X-AXIS  LIMITS  —  0.00*000.  1.00*081* 

THE  FOLLOWING  NODES  WILL  BE  PLOTTE8  -« 
10* 


The  plots  are  given  in  Fig.  B2. 


T-AXIS  LIMITS  —  70*002. -*i50*002 
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temperature  ?*ots  for  CINDA  Sample  Problem  IB.  Plots  are  shown  in  the  order  in  which  they  are  plotted.  Time  and 
temperature  values  will  be  in  the  E  format,  rather  than  in  the  F  format  as  shown. 


Fig.  B2  (Cor.t’d)— Time  vs  temperature  plots  for  CINDA  Sample 
Problem  IB.  Plots  are  shown  in  the  order  in  which  they  are  plotted. 
Time  and  temperature  values  will  be  in  the  E  format,  rather  than  in 
the  F  format  as  shown. 


